t/2 


nD-H158  170  EHPLOVNENT  OF  ADAPTIVE  LEARNING  TECHNIQUES  FOR  THE 

DISCRIHINATION  OF  ACOU.  .  <U>  GENERAL  ELECTRIC  CORPORATE 
RESEARCH  AND  DEVELOPHENT  SCHENECTA.  .  J  H  ERKES  ET  AL. 
UNCLASSIFIED  DEC  84  84SRDee2  Ne8014-82-C-2O21  F/G  2071 


NL 


EMPLOYMENT  OF  ADAPTIVE  LEARNING  TECHNIQUES 
FOR  THE  DISCRIMINATION  OF  ACOUSTIC  EMISSIONS 


rs 


< 

I 

O 

< 


J.W.  Erkes  and  K.C.  Tam 

General  Electric  Corporate  Research  and  Development 

S.R.  Mannava 

General  Electric  Turbine  Technology  Laboratory 

J.F.  MacDonald  and  H.A.  Scarton 
Rensselaer  Polytechnic  Institute 


PHASE  II  FINAL  REPORT  ^ 
Contract  N00014-82-C-2031 


December  1984 


Prepared  by 


J.W.  Erkes,  Project  Manager 
General  Electric  Company 
Corporate  Research  and  Development 
Schenectady,  New  York  12345 


Prepared  for 

H.  Chaskelis,  Project  Manager 
Naval  Research  Laboratory 
4555  Overlook  Avenue,  S.W. 
Washington,  DC  20375 


9^ 

DSTUBUnON  STATEMENT  A 

B 

Apikov*4  tot  public  icUqmi 
DiatiibvtioB  Unlimited  _ 

¥ 

84SRD002 


85  01  30  006 


UNCLASSIFIED 


SeCUniTV  CLASSIFICATION  OF  THIS  PAGE 


o-Aison^ 


REPORT  DOCUMENTATION  PAGE 


1*  REPOST  security  classification 

Unclassified 


2a.  security  classification  AUTHORITY 

N/A 


2b.  OECUASSlFICATtON/OOWNGRAOING  SCHEDULE 

N/A 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBERIS) 

84SRD002 


ID.  RESTRICTIVE  MARKINGS 

None 


S.  MONITORING  ORGANIZATION  REPORT  NUM8ER(S) 


6a.  NAME  OF  PERFORMING  QRGANIZATIC 

General  Electric  Company  / 


b.  OFFICE  SYMBOL  7a.  NAME  OF  MONITORING  ORGANIZATION 
(if  appJicadto/ 


6c.  AOORESS  (City.  St9t€  and  ZIP  Codai 

General  Electric  Company 
Corporate  Research  and  Development 
1  River  Road 
Schenectady,  NY  12345 


8a.  NAME  OF  FUNOING/'SPONSORING 

organization 

Naval  Research  Laboratory 


8c.  AOORESS  (City.  Stata  and  ZIP  Coda) 


7b.  AOORESS  (City,  Stata  and  ZIP  Coda/ 


8b.  OFFICE  symbol  9.  PROCUREMENT  INSTRUMENT  t  DENTt  IM  CATION  NUMBER 
(if  apdiicabiat 

N00014-82-C-2031  ^ 


4555  Overlook  Ave.,  S.W. 
Washington,  DC  20375 


11.  title  fineiuda  Sacunty  ciautfUattoni  Employment  of  Adaptive  Learning 
Techniques  for  the  Discrimination  of  Acoustic  Emissions  (Unclassified) 


10.  SOURCE  OF  FUNDING  NOS. 


PROGRAM 

Element  no. 


WORK  UNIT 
NO. 


12.  personal  AUTHOR(S) 


Erkes,  J.W.;  Tam,  K.C.;  Mannava,  S.R.;  MacDonald,  J.F.;  Scarton,  H.A. 


13a  type  OF  report  138.  time  covered  14.  date  of  report  (Vr.  .V/o..  Day;  15.  PAGE  COUNT 

Phase  II  Final  Report  from  Nov  82  to  Nov  83  December  1984 


COSATI  COOES 


18.  SUBJECT  TERMS  'Continua  on  rtvarse  if  neetsaary  and  identify  by  block  number/ 

'>  Acoustic  Emission,  Digital  Signal  Processing,  Adaptive  Learning, 
Homomorphic  Deconvolution, 


19.  abstract  iContinue  on  reverse  tf  necessary  and  identify  by  block  number/ 

V 

Under  Phase  I  of  this  contract,  several  new  acoustic  emission  techniques  were  developed.  These  techniques 
showed  promise  as  a  means  of  adaptively  learning  and  removing  multimode  and  multipath  effects  from  acoustic  em¬ 
ission  signals  in  a  preprocessing  step  prior  to  source  characterization.  In  this  Phase  II  effort,  software  was  developed 
to  implement  the  techniques,  and  a  series  of  experiments  was  carried  out  to  assess  the  effectiveness  and  practicality 
of  these  techniques  on  real  signals.  Analysis  of  the  data  revealed  that  reasonably  accurate  transfer  functions  could 
be  determined  adaptively,  when  a  flat  response  wide-hand  transducer  was  employed,  and  when  relatively  short 
record  lengths  were  used.  Available  commercial  transducers  with  the  requisite  frequency  response  are  quite  fragile, 
however,  and  are  far  too  insensitive  for  practical  use.  Pattern  recognition  techniques  were  also  explored  as  a  means 
of  characterizing  acoustic  emission  sources;  specifically,  details  of  the  pulse  microstructure  were  examined  for  evi¬ 
dence  of  characterizable  features.  No  such  evidence  was  found  in  either  “raw”  (reverberation  dominated)  or  pro- 
cesseddata.  A.y:.  :  - 


20.  OISTRI8UTION;AVAILABILITY  of  abstract 
UNCUASSIFIED/UNLIMITEO  X  SAME  AS  RPT.  Z  OTIC  USERS 


22a.  NAME  OF  RESPONSIBLE  iNOIVIOUAL 

H.  Chaskelis 


21.  abstract  security  CLASSIFICATION 


Unclassified 


22b  TELEPHONE  NUMBER 
•include  Area  Codei 


22c  office  symbol 


(202)  767-3613 


DD  FORM  1473,  83  APfl 


EDITION  OP  1  JAN  73  IS  OBSOLETE. 


_ Unclassified _ 

security  CLASSIFICATION  OF  THIS  PAGE 


TABLE  OF  CONTENTS 


Section  Pase 

1  INTRODUCTION  AND  SUMMARY  .  1-1 

1.1  Introduction .  1-1 

1.2  Signal  Processing .  1-1 

1.3  Sensor  Evaluation  . . .  1-2 

1.4  Experimental  Confirmation  .  1-2 

1.5  Summary .  1-2 

2  THE  THEORY  REVISITED  .  2-1 

2.1  Generalized  Time  Distribution  Function 

for  Multiple- Event  Waveforms  .  2-1 

2.2  Single-Event  Approximation  .  2-1 

2.3  Multiple  Events  .  2-2 

2.3.1  Multiple  Events  with  Occasional  Single  Events .  2-2 

2.3.2  Overlapping  Multiple  Events  .  2-2 

2.4  Exponential  Weighting .  2-4 

2.5  Bandpass  Mapping . ; .  2-7 

3  EXPERIMENTAL  APPARATUS  .  3-1 

3.1  Hardware  .  3-2 

3.1.1  Specimens .  3-2 

3.1.2  Transducers  .  3-3 

3.1.3  Analog  Signal  Conditioning  Equipment .  3-3 

3.1.4  Transient  Waveform  Recorders  .  3-3 

3.1.5  Stress-Corrosion  Event  Monitor  . .  3-4 

3.2  Software  .  3-4 

3.2.1  Data  Acquisition .  3-4 

3.2.2  Analysis .  3-5 

4  EXPERIMENTAL  RESULTS  .  4-1 

4.1  Calibration  .  4-1 

4.2  Effect  of  Record  Length  .  4-4 

4.3  Dependence  on  the  Width  of  the  Cepstral  Filter  .  4-6 

4.4  Dependence  on  Alpha .  4-6 

4.5  Fourier  Deconvolution  .  4-9 

4.6  Effects  of  a  Limited-Frequency  Band .  4-14 

4.7  Pattern  Recognition .  4-16 

5  CONCLUSION  .  5-1 

6  REFERENCES  .  6-1 


TABLE  OF  CONTENTS  (Confd) 

Section 

APPENDIX  A  -  BIOMATION  8100  DRll-C  INTERFACE 
APPENDIX  B  -  ILS  SOFTWARE  MODULES 

APPENDDC  C  -  PATTERN  RECOGNITION  SCATTER  PLOT  PRODUCTION 


DTIC 


LIST  OF  ILLUSTRATIONS 


■ 

* 

m 

i-'- 


t 

I 

* 

).  •• 


t*  'j 

i'- 


t 


Figure  Page 

1  The  delta-shaped  time  distribution  function  obtained  by  averaging  the  time 
distribution  function  of  a  large  number  of  multiple-event  waveforms  with  one 

of  the  events  in  each  waveform  aligned  .  2-4 

2  The  relation  between  record  length  and  frequency  sampling  interval  .  2-5 

3  Removing  zeroes  in  z-transform  by  exponential  weighting .  2-6 

4  Modified  homomorphic  characteristic  system  including  exponential  weighting  .  2-6 

5  Bandpass  mapping  operation  .  2-7 

6  Characteristic  system  for  bandpass  homomorphic  systems  .  2-8 

7  Block  diagram  of  the  data  acquistion  system .  3-1 

8  The  experimental  setup,  showing  the  thick  test  plate,  the  Biomation  transient 

recorders,  the  analog  instrumentation,  and  VAX  11/750  control/analysis  computer  ....  3-2 

9  A  typical  lead-break  waveform  recorded  by  a  NBS  conical  transducer .  4-2 

10  Another  lead-break  waveform  recorded  by  a  NBS  conical  transducer  located 

54  cm  from  the  one  used  in  Figure  9  .  4-2 

11  The  cepstrum  of  the  waveform  in  Figure  10 . 4-3 

12  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  10 .  4-4 

13  The  improved  result  obtained  through  the  use  of  exponential  weighting;  a  -  0.9955  ..  4-5 

14  The  theoretical  seismic  surface  pulse  according  to  Pekaris .  4-6 

15  A  lead-break  waveform  recorded  by  a  NBS  conical  transducer  located  45  cm 

from  the  source  .  4-7 

16  A  lead-break  waveform  recorded  by  a  NBS  conical  transducer  located  45  cm 

from  the  source .  4-7 

17  A  lead-break  waveform  recorded  by  a  NBS  conical  transducer  located  30  cm 

from  the  source .  4-8 

18  A  lead-break  waveform  recorded  by  a  NBS  conical  transducer  located  15  cm 

from  the  source .  4-8 

19  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  15, 

using  all  of  the  8096  samples  as  input;  a  -  0.9985  .  4-9 

20  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  16, 

using  all  of  the  8096  samples  as  input;  a  -  0.9985  .  4-10 

21  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  17, 

using  all  of  the  8096  samples  as  input;  a  -  0.9985  .  4-10 

22  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  18, 

using  all  of  the  8096  samples  as  input;  a  -  0.9985  .  4-11 


iK 


V 


LIST  OF  ILLUSTRATIONS  (Cont’d) 


Figure  Page 

23  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  15, 

using  only  the  first  2048  samples  as  input;  a  -  0.9985  .  4-12 

24  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  16, 

using  only  the  first  2048  samples  as  input;  a  -  0.9985  .  4-13 

25  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  17, 

using  only  the  first  2048  samples  as  input;  a  -  0.9985  .  4-13 

26  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  18, 

using  only  the  first  2048  samples  as  input;  a  -  0.9985  .  4-14 

27  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  10, 

a  -  0.999,  width  of  cepstral  filter  -  ±2.5  fis  .  4-15 

28  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  10, 

a  -  0.999,  width  of  cepstral  filter  -  ±4.5  .  4-16 

29  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  10, 

a  -  0.999,  width  of  cepstral  filter  -  ±9.5  fis  .  4-17 

30  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  15, 

a  -  0.9985,  width  of  cepstral  filter  -  ±2.5  .  4-18 

31  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  15, 

a  -  0.9985,  width  of  cepstral  filter  -  ±4.5  .  4-18 

32  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  15, 

a  -  0.9985,  width  of  cepstral  filter  -  ±9.5  /its  .  4-19 

33  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  15, 

a  -  0.9985  .  4-19 

34  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  15, 

a  -  0.9980  .  4-20 

35  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform  in  Figure  15, 

a  -  0.9975  .  4-20 

36  A  typical  large-angle  Pentel  lead-break  waveform;  the  angle  between  the  lead 

and  the  surface  is  65® .  4-21 

37  A  typical  small-angle  Pentel  lead-break  waveform;  the  angle  between  the  lead 

and  the  surface  is  40® .  4-21 

38  A  Pentel  lead-break  waveform  with  secondary  Pentel  tip  impact  at  a  large  angle  .  4-22 

39  A  200-point  sample  of  a  fixture-generated  lead-break  waveform 

near  the  initial  impulse  .  4-22 

40  A  200-point  sample  of  another  fixture-generated  lead-break  waveform 

near  the  initial  impulse  .  4-23 


LIST  OF  ILLUSTRATIONS  (Cont’d) 

Figure  Page 

41  Another  200-point  sample  of  the  waveform  in  Figure  39  near  the  later  part 

of  the  waveform .  4-23 

42  Another  200-point  sample  of  the  waveform  in  Figure  40  near  the  later  part 

of  the  waveform .  4-24 

43  The  impulse  response  obtained  by  high-time  pass  filtering  the  cepstrum  of  the  average 

of  33  mechanically  generated  lead-break  waveforms;  a  -  0.999  .  4-25 

44  One  of  the  33  mechanically  generated  lead-break  waveforms  used  in  the  averaging  ....  4-26 

45  The  result  of  deconvolving  the  waveform  in  Figure  44  by  the  impulse  response 

in  Figure  43  .  4-27 

46  An  elliptical  frequency  filter  used  to  reduce  the  high-frequency  noise  in  the  result 

in  Figure  45;  -60  dB  cut-off  at  700  kHz .  4-28 

47  The  result  of  filtering  the  waveform  in  Figure  45  with  the  elliptical  filter  in  Figure  46  4-29 

48  Relative  amplitude  of  the  stress  corrosion  events  and  the  lead-break  events .  4-29 

49  Power  spectrum  of  a  stress-corrosion  event  captured  by  a  resonance  transducer  .  4-30 

50  The  average  of  31  stress-corrosion  waveforms  with  their  initial  impulses  lined  up 

at  the  same  location;  the  waveforms  were  captured  with  a  resonance  transducer  .  4-31 

51  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  averaged  waveform 

in  Figure  50;  no  band-pass  mapping  .  4-32 

52  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  averaged 

waveform  in  Figure  50;  band-pass  mapping  between  100  and  300  kHz .  4-33 

53  Scatter  plots  of  the  first  two  principal  components  of  the  frequency  components 

of  two  sets  of  events  after  low-time  pass  filtering  in  the  cepstral  domain .  4-34 

54  Scatter  plots  of  the  first  two  principal  components  of  the  frequency  components 

of  the  two  sets  of  events  in  Figure  53  without  filtering  in  the  cepstral  domain .  4-35 

55  Scatter  plots  of  the  first  two  principal  components  of  the  frequency  components 

of  three  sets  of  events  after  low-time  pass  filtering  in  the  cepstral  domain  .  4-36 

56  Scatter  plots  of  the  first  two  principal  components  of  the  frequency  components 

of  the  three  sets  of  events  in  Figure  55  without  filtering  in  the  cepstral  domain  .  4-37 


Section  1 


INTRODUCTION  AND  SUMMARY 


1.1  Introduction 

Until  recently,  workers  in  the  acoustic  emission  (AE)  research  community  appear 
to  have  been  unaware  of  a  number  of  techniques  applied  successfully  in  other  research 
areas  (including  sonar,  seismology,  and  speech  processing)  where  certain  advanced  sig¬ 
nal  processing  techniques  are  frequently  used  to  deal  with  distortion  effects  imposed 
on  the  signals  of  interest.  For  example,  volume  reverberation  and  multimode  propaga¬ 
tion  effects  are  common  problems  in  these  applications,  and  much  effort  has  been  ex¬ 
pended  in  the  development  of  advanced  techniques  to  compensate  for  them.  Since 
these  effects  also  severely  distort  real  AE  signals,  those  advanced  techniques  are  ex¬ 
pected  to  be  useful  in  improving  the  performance  and  reliability  of  AE-based  nonde¬ 
structive  evaluation  (NDE)  systems.  This  report  describes  the  results  of  a  research 
program  that  seeks  to  capitalize  on  those  successful  research  efforts,  extending  them 
and  applying  them  to  the  problems  of  the  detection,  location,  and  characterization  of 
AE  events. 

1.2  Signal  Processing 

Current  AE  analysis  techniques  make  limited  use  of  the  phase  information  in  the 
received  signal  and  rely  on  a  variety  of  incoherent  techniques  for  source  identification 
and  location;  pulse  energies,  rise  times,  ringdown  times,  and  amplitude  distributions 
are  commonly  used  signal  parameters.  Although  these  parameters  often  yield  useful 
information  about  the  nature  and  location  of  the  AE  source,  they  do  have  serious  lim¬ 
itations.  In  particular,  they  tend  to  be  less  useful  when  the  signal-to-noise  (S/N)  ratio 
is  poor,  when  complex  structural  features  produce  confusing  reflections,  or  when  mul¬ 
tiple  sources  are  present.  Unfortunately,  the  most  important  application  areas  for 
acoustic  emissions  often  suffer  badly  from  just  these  problems.  Specifically,  on-line 
AE-based  NDE  systems  often  must  deal  simultaneously  with  adverse  S/N  ratios,  com¬ 
plex  structures  producing  serious  multipath/muitimode  interference,  as  well  as  un¬ 
known  (often  multiple)  source  locations. 

Coherent  methods  (which  make  use  of  the  phase  information)  provide  some  basis 
for  hope  in  dealing  with  these  complex  problems.  Unfortunately,  linear  methods  to 
compensate  for  these  effects  have  not  proven  to  be  particularly  useful.  In  light  of  the 
role  of  mode  conversion  and  multimode  propagation  in  signal  distortion,  this  is 
perhaps  not  surprising.  Nonlinear  methods  offer  considerably  more  promise  and,  in 
fact,  have  been  used  with  much  success  in  dealing  with  similar  problems  in  other 
fields.  Homomorphic  deconvolution  especially  has  been  used  with  considerable  suc¬ 
cess  to  eliminate  multipath  and  multimode  distortion  in  seismic,  speech,  and  audio 
processing  applications  where  S/N  ratios  are  relatively  good.  A  series  of  potentially 
useful  homomorphic  signal  processing  approaches  were  developed  under  Phase  1  of 
this  contract  and  are  described  in  detail  in  the  Phase  I  final  report.  This  report  will 
focus  on  the  experimental  evaluation  of  the  potential  use  these  techniques  may  ulti¬ 
mately  find  in  practical  acoustic  emission  problems  and,  in  particular,  evaluate  their 


potential  for  compensating  for  transducer  and  structural  resonances.  Reliable  mul¬ 
timode  and  multipath  compensation,  if  achievable,  can  provide  the  technical  basis  for 
automated  on-line  AE  monitoring  of  complex  structures  for  cracks;  if  the  structural 
complexities  can  be  reduced  or  eliminated,  it  may  be  practical  to  monitor  crack  growth 
rates,  to  carry  out  crack  severity  assessments,  and  to  locate  cracks. 

In  other  situations,  particularly  when  the  S/N  ratio  is  poor,  homomorphic 
methods,  especially  those  based  on  cepstral  techniques,  are  less  useful  and  suffer  from 
problems  with  accurate  phase  unwrapping.  Under  these  conditions,  adaptive  nonlinear 
methods  may  be  very  applicable.  The  experimental  work  reported  here  focused  on  the 
deconvolution  problem  for  situations  where  a  relatively  good  S/N  ratio  was  available. 

1.3  Sensor  Evaluation 

The  most  commonly  used  AE  transducers — resonant  piezoelectric  transducers — 
seriously  distort  the  microstructure  of  AE  signals  through  transducer  ringdown. 
Although  the  nonlinear  signal  processing  approaches  outlined  above  will  to  a  large  ex¬ 
tent  compensate  for  those  problems,  a  better  solution  would  be  to  use  a  transducer 
without  massive  intrinsic  phase  distortion.  In  this  project  we  studied  the  performances 
of  both  the  resonant  transducers  and  the  conical  transducers.  The  conical  transducers 
used  were  commercially  available  piezoelectric  sensors  based  on  ideas  developed  and 
demonstrated  by  Eitzen  et  al.'  at  the  National  Bureau  of  Standards.  They  have  rela¬ 
tively  flat  frequency  response  up  to  1  MHz.  Such  improved  frequency  response  is 
achieved  by  reducing  the  contact  area  of  the  active  element,  by  creating  a  smooth 
acoustical  impedance  transition  to  the  backing  material,  and  by  increasing  the  backing 
material  volume.  Because  of  the  reduced  contact  area,  however,  the  detected  signals 
are  much  reduced  in  amplitude  compared  to  the  resonant  transducers,  and  very  poor 
S/N  ratios  can  be  a  serious  problem  in  practical  experimental  situations. 

1.4  Experimental  Confirmation 

Finally,  as  the  major  element  of  the  Phase  II  contract  effort,  a  series  of  experi¬ 
ments  of  gradually  increasing  complexity  were  planned  and  performed  to  confirm  and 
evaluate  the  effectiveness  of  the  signal  processing  algorithms  and  techniques  devel¬ 
oped.  The  experiments  were  carried  out  by  General  Electric  at  the  Research  and 
Development  Center  and  at  the  Materials  and  Processes  Laboratory,  in  Schenectady, 
N.Y.,  under  the  joint  direction  of  the  co-investigators. 

1.5  Summary 

This  project  represents  a  radically  new  approach  to  the  acquisition  and  analysis  of 
AE  signals.  The  approach  is  based  on  the  use  of  new  signal  processing  methods  and 
sensors  and  focuses  on  the  removal  of  multipath  and  multimode  distortion  from  the 
signals. 


Section  2 


THE  THEORY  REVISITED 


The  analytical  basis  for  this  work  is  described  in  detail  in  the  Phase  I  final  report, 
and  will  not  be  repeated  here.  It  should  be  noted,  however,  that  in  the  course  of  car¬ 
rying  out  the  experimental  portion  of  Phase  II,  it  was  discovered  that  some  of  the 
theoretical  results  developed  in  Phase  I  of  this  project  needed  to  be  extended,  and,  in 
some  instances,  modified.  Some  of  these  changes  came  about  in  the  course  of  analyz¬ 
ing  the  experimental  data  and  in  learning  more  about  actual  experimental  uncertain¬ 
ties.  Other  changes  arose  out  of  a  more  careful  examination  of  the  physical  assump¬ 
tions  behind  the  application  of  these  techniques  to  the  AE  problem.  The  topics 
modified  or  added  include  exponential  weighting,  bandpass  mapping,  and  the  proper 
choice  of  the  time  distribution  function  in  analyzing  the  multipie-event  waveforms. 
Detailed  descriptions  of  these  changes,  and  the  reasons  behind  them,  are  outlined 
below. 


2.1  Generalized  Time  Distribution  Function  for  Multiple- Event  Waveforms 

As  was  shown  in  the  Phase  I  report,  the  measured  signals  y(t)  in  acoustic  emission 
experiments  can  be  written  in  the  form 


y(t) 


T  j) 


hfit) 


g{t)+  nit) 


(1) 


where 


•  =  the  operation  of  convolution 

h,U)  =  the  ringdown  impulse  response  due  to  the  boundary  reverberation 
and  the  transducer  response 

hpiU—Tj)  =  the  AE  impulse  pulse  emitted  at  the  source  at  time  t, 

git)  =  the  measurement  time  gate 
nit)  =  the  background 

The  duration  of  the  measurement  time  gate  git)  is  assumed  to  be  much  longer  than 
ihose  of  hp,  and  h,  and  therefore  will  be  omitted  from  the  rest  of  this  report.  This 
point  will  be  treated  in  more  detail  is  Section  4.2. 

2.2  Single-Event  Approximation 

Single  events  represent  a  much  simpler  case,  and  in  the  case  of  a  solitary  event. 
Equation  1  reduces  to 

yit)  -  hpjit)  *  hrit)+  nit)  (2) 

if  the  noise  component  nit)  is  negligible.  Equation  2  can  be  treated  by  the  conven¬ 
tional  homomorphic  deconvolution  technique  to  recover  the  signal  /tp,(r).  Under  this 
approach,  the  measurement  yit)  is  transformed  to  cepstral  domain  by  the 
homomorphic  operation  IFT(LOG(FT0’(r)))).  This  composite  operation  converts  the 
convolution  operation  into  the  addition  operation.  The  cepstrum  of  yit)  is  the  sum  of 
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those  of  hpjU)  and  hfU),  and  hence  the  two  cepstra  can  be  separated  by  linear  filter¬ 
ing. 

This  technique  will  work  especially  well  when  the  cepstra  of  hj„U)  and  h,(t)  are 
well  separated  in  cepstral  period.  In  A£  measurements,  these  conditions  typically  are 
met;  typical  AE  /?p,s  are  narrow  impulses  with  durations  on  the  order  of  a  few  mi¬ 
croseconds,  while  the  h,%  have  substantial  low-frequency  energy  components  and  are 
much  longer  in  duration.  Since  the  cepstrum  of  a  delta-like  function  is  also  a  delta-like 
function,  the  cepstra  of  the  /ip,s  are  concentrated  near  the  origin  in  the  cepstral 
domain.^-^  In  general  the  cepstrum  of  A,  is  more  spread  out;  so  the  criterion  of  separa¬ 
bility  of  cepstra  is  approximately  satisfied.  One  way  to  further  reduce  the  contribution 
of  the  cepstrum  of  h,  near  the  origin  is  to  convert  into  a  minimum-phase  sequence 
by  using  the  technique  of  exponential  weighting,  which  will  be  described  in  detail  in 
Section  2.4. 

Once  the  cepstra  of  and  KU)  are  separated,  they  can  be  operated  on  by  the 
inverse  homomorphic  transform  to  yield  either  hpjU)  and  respectively.  The 

recovered  pulse  shape  can  be  analyzed  for  features  associated  with  the  particular 
class  of  events  characterizing  the  source. 

2.3  Multiple  Events 

In  general,  AE  events  do  not  occur  singly,  but  rather  in  associated  bursts;  in  some 
structures,  the  events  may  be  easily  visible  as  single  pulses,  but  in  large  reverberant 
structures,  the  individual  pulses,  extended  by  reverberation,  will  overlap  and  be 
difficult  to  distinguish  from  one  another. 

2.3.1  Multiple  Events  with  Occasional  Single  Events 

On  the  other  hand,  provided  that  the  AE  events  are  associated  with  a  localized 
crack  formation  and  are  statistically  stationary,  and  further  provided  that  a  single  re¬ 
verberated  AE  waveform  is  available  for  analysis,  the  derived  ringdown  function  h^{t) 
can  be  used  to  dereverberate  the  multiple-event  waveforms  (if  they  occur)  from  the 
same  source,  through  simple  Fourier  inversion.  These  dereverberated  waveforms 
would  reveal  the  original  impulse  time  sequence  of  the  overlapping  individual  events, 
and  hence  similar  data  from  another  transducer  located  at  a  different  spatial  location 
could  be  analyzed  jointly  via  cross-correlation  techniques  to  yield  information  on  the 
time  of  arrival  and  hence  location  of  the  events.  Similarly,  the  rate  of  occurrence  of 
these  waveforms  could  be  easily  determined  and  could  be  used  in  estimating  the  ac¬ 
tivity  strength  and  otherwise  characterizing  the  source. 

2.3.2  Overlapping  Multiple  Events 

More  generally,  in  dealing  with  multiple-event  waveforms,  one  has  to  resort  to 
Equation  1.  Compared  to  Equation  2,  Equation  1  has  two  complications:  the  constancy 
of  the  impulse  shape  hpj  and  the  time  distribution  of  the  events.  Assuming  the 
different  events  all  have  the  common  pulse  shape  hp{t)  and  only  differ  in  their  ampli¬ 
tudes  n,,  then  we  have 

y{t)  -  hp(t)  *  pit)  *  hrit)+n(t)  (3) 

where  hpit)  is  the  AE  pulse  emitted  at  time  r,  and  pit)  -  ^a(/)8(t-T,)  is  the  time 


distribution  of  the  events  weighted  with  amplitudes  corresponding  to  those  of  the 
different  events. 

The  cepstrum  of  the  measured  signal  >»(/)  is  then  the  sum  of  the  cepstra  of 
pit),  and  hrit).  If  the  pulse  train  pit)  has  been  converted  to  minimum  phase,  it  can 
been  shown  that  its  cepstrum  is  nonzero  only  for  positive  times  greater  than  or  equal 
to  the  interval  between  the  first  two  arrivals,^  and  thus  is  well  separated  from  that  of 
hpit)  if  the  intervals  between  the  impulses  are  not  too  close  to  each  other.  In  this 
case,  hpit)  can  again  be  recovered  by  homomorphic  deconvolution. 

(If  the  different  events  do  not  have  the  same  pulse  shape,  the  signals 
'J^Ojhpiit  —  Ti)  contained  in  the  measurement  yit)  cannot  be  simplified  to  a  convolu- 
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tion  of  the  pulse  train  pit)  with  a  common  pulse  shape  hpit),  as  in  Equation  3.  Since 
it  is  difficult  to  estimate  the  cepstra  of  such  complex  signals,  cepstral  analysis  is  not 
practical.) 

On  the  other  hand,  since  the  cepstra  of  h,it)  and  pit)  usually  overlap,  h,if)  can¬ 
not  be  separated  out.  (Note  that  this  result  is  distinctly  different  from  the  result  re¬ 
ported  in  the  Phase  I  study,  where  the  assumption  that  an  exponential  distribution 
function  holds  for  an  ensemble  of  overlapping  AE  pulses  lead  to  a  contrary  conclu¬ 
sion.  A  more  careful  analysis  of  the  physical  situation  suggests  that  an  exponential  dis¬ 
tribution  is  not  appropriate  for  an  overlapping  ensemble  of  AE  pulses.)  Fortunately 
there  is  another  way  out  of  this  dilemma.  The  method  works  on  the  fact  that  if  pit)  is 
in  fact  a  delta-like  function,  then  the  cepstrum  of  pit)  is  also  a  delta-like  function  and 
thus  is  well  separated  from  that  of  /?,(/).  One  way  to  simulate  this  situation  is  by 
averaging  a  large  number  of  waveforms  with  one  of  the  events  in  each  waveform 
aligned  at  about  the  same  location.  A  simple  and  practical  way  of  accomplishing  this, 
for  example,  is  to  use  a  transient  waveform  capturing  device  set  to  trigger  on  a  strong 
impulse.  This  will  result  in  a  series  of  captured  waveforms  where  the  first  strong  event 
in  each  waveform  which  exceeds  the  triggering  level  always  occurs  at  about  the  same 
location  /q  in  the  waveform  records.  Under  these  circumstances,  the  time  distribution 
function  can  be  approximated  by 

pit)  ~  Sit- to)  + b  (4) 

where  b  represents  the  averaged  time  distribution  of  the  other  events  in  the 
waveforms. 

The  procedure  is  illustrated  in  Figure  1.  If  the  number  of  waveforms  N  used  in  the 
averaging  is  large  enough,  b  will  approach  a  uniformly  flat  distribution.  The  magnitude 
of  b  decreases  with  the  number  N ,  approaching  zero  as  N  tends  to  infinity.  In  this 
case  pit)  is  basically  a  delta  function  with  a  cepstrum  concentrated  at  the  origin,  and 
thus  will  not  interfere  with  the  cepstrum  of  h^it). 

Note  that  if  there  were  no  aligning  in  the  averaging,  all  the  impulses  would  be  dis¬ 
tributed  randomly  in  location,  and  pit)  would  approach  a  uniformly  flat  distribution. 
Since  the  convolution  of  a  flat  distribution  with  any  function  is  also  a  flat  distribution, 
the  averaged  waveform  yit)  is  also  a  flat  distribution  and  therefore  contains  no  infor¬ 
mation  at  all. 
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Figure  1.  The  delta-shaped  time  distribution  function  obtained  by  averaging  the 
time  distribution  function  of  a  iarge  number  of  multiple-event  waveforms 
with  one  of  the  events  in  each  waveform  aligned 

2.4  Exponential  Weighting 

One  of  the  steps  in  cepstral  analysis  is  phase  unwrapping,  which  is  the  process  of 
making  the  phase  0  of  the  Fourier  transform,  rexp((«),  of  the  time  series  waveform 
{a(n)}  into  a  continuous  function.’  However,  if  some  Fourier  components  are  zero, 
their  phase  9  would  be  undefined,  and  phase  unwrapping  would  fail.  This  problem  is 
especially  serious  with  waveforms  of  long  record  length.  The  reason  for  this  is  that 
long  record  length  corresponds  to  fine  sampling  in  the  Fourier  domain,'^  and  hence 
waveforms  with  long  record  length  are  more  likely  to  pick  up  zero  Fourier  com¬ 
ponents  than  those  with  short  record  length.  This  phenomenon  is  illustrated  in  Fig¬ 
ure  2. 

Now  the  Fourier  transform  of  a  waveform  corresponds  to  the  r-transform  of  the 
function  at  the  unit  circle  |r|=  1,  where  the  r-transform  A(z)  of  a  sequence  a(n)  is 
defined  as: 

(z)  -  X  o{n)2^". 


where  z  is  a  complex  variable. 
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Figure  2.  The  relation  between  record  length  and  frequency  sampling  interval 

One  way  to  get  around  the  phase  wrapping  problem  when  some  Fourier  com¬ 
ponents  are  zero  is  by  using  exponential  weighting,’  The  original  waveform  |a(«)}  is 
weighted  with  the  series  (a"}  before  z-transformation: 

where  a  is  a  real  number  less  than  1.  Exponentially  weighting  a  function  by  |a'')  has 
the  effect  of  scaling  the  magnitude  of  the  zeroes  of  the  z-transform  of  the  function  by 
a,  as  illustrated  in  Figure  3.  Thus  the  zeroes  in  the  Fourier  transform  of  a  function 
can  be  removed  by  exponentially  weighting  the  function  before  inputting  to  the 
homomorphic  system,  and  phase  unwrapping  can  be  successfully  carried  out. 

The  use  of  exponential  weighting  can  also  improve  the  separability  of  the  cepstra. 
If  a  is  small  enough  so  that  all  the  zeroes  of  the  reverberation  sequence  are  scaled  to 
lie  within  the  unit  circle,  the  sequence  is  converted  to  a  minimum-phase  sequence. 
Under  this  condition  the  reverberation  sequence  will  contribute  to  the  cepstrum  only 
for  positive  times  greater  than  or  equal  to  the  interval  between  the  first  two  arrivals  in 
the  sequence,^  and  hence  the  separability  of  the  cepstra  will  be  improved.  In  contrast, 
the  structure  of  the  cepstrum  of  a  mixed-phase  sequence  is  complicated,  difficult  to 
predict,  and  frequently  unstable. 

Note  that  exponential  weighting  is  different  from  the  usual  windowing  procedure 
used  in  Fourier  transformation,  which  is  an  approximation  as  far  as  convolution  is 
concerned.  In  general,  convolution  is  not  conserved  under  multiplication  by  a  win¬ 
dowing  function:  the  product  of  a  windowing  function,  win),  with  the  convolution  of 


|Z|  =  1,  UNIT  CIRCLE* 
FOURIER  TRANSFORMATION 


X,,  xj  ZEROES  OF  z  TRANSFORM  OF  x(n) 
y,,  y2  ZEROES  OF  z  TRANSFORM  OF  o"  x(n) 

IVii  =  o  X, 

'  Vzi  =  o  X2 

Figure  3.  Removing  zeroes  in  z-tronsform  by  exponentiai  weighting 

two  functions,  /(«)  and  gift),  is  different  from  the  convolution  of  w{n)f{n)  and 
w(n)g(«)'. 

tv(rt)|/(rt)  *  ^  |w(n)/(/i)j  •  |»v(/i)^(/i)| 

Convolution  will  be  conserved,  however,  if  w(n)  is  in  the  form  a": 

a'-j/f/i)  *  ^(«))  -  (a”/(«))  •  («"?(«)) 

The  output  from  the  inverse  homomorphic  system  is  deweighted  by  the  sequence 
{l/a").  The  modified  homomorphic  characteristic  system  including  exponential  weight¬ 
ing  is  illustrated  in  Figure  4. 


2.5  Bandpass  Mapping 

If  the  signals  one  wishes  to  analyze  have  intrinsic  bandpass  characteristics,  as,  for 
example,  the  data  taken  with  resonant  transducers,  homomorphic  analysis  cannot  be 
accomplished  by  means  of  full-band  homo.morphic  systems.  In  fact,  the  analysis  can¬ 
not  be  performed  on  the  unit  circle  of  the  z-plane,  since  then  the  logarithm  would  be¬ 
come  unbounded  in  the  frequency  bands  with  zero  energy.  Neither  can  it  be  per¬ 
formed  off  the  unit  circle,  since  the  z-transform  of  such  signals  does  not  converge 
anywhere  on  the  z-plane,  but  on  the  unit  circle.  Thus  the  exponential  weighting  pro¬ 
cedure  used  to  remove  the  zeroes  of  a  signal  off  the  unit  circle  cannot  be  applied  to 
the  class  of  bandpass  signals.  Since  real  data  always  contain  some  out-of-band  noise, 
and  can  only  have  a  countable  number  of  zeroes  on  the  unit  circle,  sometimes  it  ap¬ 
pears  possible  to  employ  the  exponential  weighting  procedure  to  remove  such  zeroes 
from  the  unit  circle.  Such  an  approach,  however,  is  inherently  ill-conditioned,  and 
often  leads  to  erroneous  results. 

One  way  to  solve  this  problem  involves  a  restriction  on  the  domain  of  the  logarith¬ 
mic  mapping  to  encompass  only  the  passband  of  the  input.  This  approach  may  be  con¬ 
veniently  formulated  in  terms  of  a  frequency  scaling  operation  that  shifts  and  stretches 
the  signal’s  passband  to  occupy  the  entire  frequency  domain,  as  illustrated  in  Figure  S. 
The  result  of  this  operation  is  then  a  full-band  sequence,  which  may  be  analyzed  using 
full-band  homomorphic  systems.  The  analyzed  output  from  the  homomorphic  system 
undergoes  the  inverse  bandpass  mapping.  The  characteristic  system  for  bandpass 
homomorphic  systems  is  illustrated  in  Figure  6. 


/  \ 
/  \ 
/  SHIFTING  &  INTERPOLATION  \ 


/  \ 


Figure  5.  Bandpass  mapping  operation 


Figure  6.  Characteristic  system  for  bandpass  homomorphic  systems 


Bear  in  mind  that  bandpass  mapping  is  not  a  panacea,  but  rather  a  method  to  re¬ 
move  the  effects  of  the  out-of-band  noise  in  phase  unwrapping.  It  does  not  add  any  in¬ 
formation  in  the  frequency  bands  with  zero  energy  and  may  not  produce  particularly 
satisfying  results,  especially  if  the  system  under  analysis  is  sharply  resonant,  with  a  rel¬ 
atively  small  range  of  useful,  information-carrying  frequencies. 


Section  3 


EXPERIMENTAL  APPARATUS 


A  VAX-based  experimental  system  was  assembled  to  provide  a  suitable  testbed  for 
the  acoustic  emission  experiments  carried  out  under  this  contract.  The  system,  shown 
in  schematic  form  in  Figure  7,  provided  a  great  deal  of  flexibility  in  the  acquisition, 
validation,  and  analysis  of  the  acoustic  emission  data.  Data  acquisition,  for  example, 
was  carried  out  under  VAX  control;  this  approach  gave  the  experimenters  access  to 
VAX-based  tools  to  monitor  the  newly  acquired  data  and  verify  its  integrity  during  the 
course  of  the  experiments.  In  addition,  the  VAX  system  incorporated  an  interactive 
digital  signal  processing  software  system,  which  was  used  extensively  in  the  analysis  of 
the  data;  the  system  was  useful  not  only  in  the  exploratory  stages  while  appropriate  al¬ 
gorithms  were  being  developed,  but  also  during  the  production  analysis  phase  when 
validated  experimental  data  were  being  batch  processed.  Subsequent  sections  of  this 
report  will  provide  details  on  the  hardware  and  software  used  in  these  experiments. 


FIGURE  8  -  BLOCK  DIAGRAM  OF  THE  DATA  ACQUISITION  SYSTEM 


Figure  7.  Block  diagram  of  the  data  acquisistion  system 


3.1  Hardware 


The  hardware  used  in  the  acoustic  emission  experiments  will  be  described  in  the 
subsequent  sections  essentially  in  the  same  order  as  the  signal  flow  through  the  sys¬ 
tem.  Figure  8  shows  the  experimental  setup;  the  thick  test  plate  can  clearly  be  seen,  as 
can  the  transducers,  the  Biomation  transient  recorders,  the  analog  instrumentation, 
and  the  VAX  11/750  control/analysis  computer  in  the  background. 


Figure  8.  The  experimental  setup,  showing  the  thick  test  plate,  the  Biomation 
transient  recorders,  the  analog  instrumentation,  and  VAX  11/750 
control/anaiysis  computer 


3.1.1  Specimens 

A  series  of  experimental  specimens  was  prepared  to  provide  a  reasonable  range  of 
expected  reverberation  conditions;  a  key  element  in  the  design  was  the  notion  that  it 
should  be  possible  to  increase  gradually,  and  in  a  controlled  manner,  the  physical  com¬ 
plexity  of  the  experiments,  so  that  new  conditions  were  added  one  at  a  time. 

3. 1.1.1  Workpieces.  Two  aluminum  plates  were  selected  to  conduct  the  experiments. 
The  first  one  was  a  relatively  thick  plate  with  dimensions  of  27  x  5  x  2  in.  The  thick 
plate  was  chosen  to  allow  clean  near-held  measurements  without  the  confusing  effects 
of  reverberations  from  the  back  wail  and  other  reflections.  It  functions  in  a  near-field 
sense  as  an  infinite  half-sphere,  and  has  known  solutions  for  the  expected  AE  due  to  a 
Pentel  pencil  lead  break.  The  first  plate  is  set  up  on  a  2  in.  thick  foam  pad  to  isolate 
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from  ambient  vibrations.  The  second  plate  had  dimensions  of  21  y.  IS  x.  1/4  in.  and 
was  chosen  to  allow  the  investigation  of  the  effects  of  plate  bending  vibrations  and 
near-field  resonance,  as  well  as  edge  reflections.  Only  the  thick  plate  was  used  for 
most  of  the  experiments  described  in  this  report. 

3. 1.1.2  C-Block  Stress  Corrosion  AE  Sources.  Notched  aluminum  C-blocks,  under 
tensile  stress  and  subject  to  a  corrosive  solution,  were  used  to  generate  natural  acous¬ 
tic  emissions.  The  notched  C-blocks  were  clamped  to  the  test  plate  and  coupled  with 
glycerine,  so  that  good  transmission  of  the  generated  AE  pulses  was  achieved.  A  sodi¬ 
um  chloride  solution  was  subsequently  used  to  induce  stress  corrosion  cracking  in  the 
notched  block.  These  sources  provided  reliable,  natural  AE  signals,  and  could  be  easily 
controlled  by  small  adjustments  in  the  stress  and/or  salt  solution  so  as  to  provide  a 
wide  range  of  natural  AE  signals,  ranging  from  weak,  individual  pulses  to  nearly  con¬ 
tinuous  trains  of  strong  pulses. 


3.1.2  Transducers 

Both  conventional  resonant  transducers  (Physical  Acoustics  Corporation,  Model 
PAC  RIS)  as  well  as  conical  low-resonance  wide-band  transducers  (Industrial  Quality 
Inc.,  Gaithersburg,  Md.)  were  used  to  detect  the  acoustic  events.  The  resonant  trans¬ 
ducers  have  a  band  frequency  response  between  100  and  300  kHz.  The  conical  trans¬ 
ducers  have  a  flat  frequency  response  up  to  1  MHz. 

3.1.3  Analog  Signal  Conditioning  Equipment 

The  signals  from  the  transducers  were  first  amplified  using  a  Tektronix  Model  502 
amplifier  with  internal  low-pass  filter  set  at  1  MHz,  or  lower,  depending  on  the  nature 
of  the  sensing  transducer.  These  signals  were  further  processed  using  a  Kronhite  pro¬ 
grammable  filter  (48  dB/octave/section)  set  up  to  low  pass  below  1  MHz.  The  output 
of  the  filter  was  then  connected  to  the  Biomation  transient  recorder  system.  A  Nicolet 
digital  Explorer  scope  was  used  for  a  visual  display  of  the  signals  prior  to  transmitting 
them  into  the  computer. 

3.1.4  Transient  Waveform  Recorders 

Biomation  8100  digital  transient  recorders  were  used  to  record  and  sample  the 
acoustic  emission  signals.  Using  a  single  Biomation  8100,  2048  samples  can  be  ob¬ 
tained  in  single-channel  mode;  in  dual-channel  mode,  two  parallel  channels  of  1024 
samples  can  be  obtained.  The  recorder  operates  at  up  to  a  maximum  sampling  rate  of 
100  MHz,  if  desired.  For  single-channel  operation,  a  2  MHz  sampling  rate  results  in 
the  acquisition  of  nearly  1  ms  of  data.  Longer  records  needed  to  verify  the  Phase  1  re¬ 
sults  were  obtained  by  ganging  together  four  of  these  systems  in  series.  Using  this  ar¬ 
rangement,  nearly  4  ms  of  continuous  data  could  be  acquired  from  a  single  AE  event, 
if  required.  A  VAX-Biomation  interface  board  was  designed  and  fabricated  to  provide 
both  automated  experiment  control  by  the  VAX,  as  well  as  a  data  path  for  the 
transmission  of  transducer  signals  from  the  Biomation  recorder  into  the  DEC  VAX 
11/750  computer  for  later  postprocessing.  The  hardware  and  software  for  this  data  ac¬ 
quisition  system  were  designed  and  put  together  during  the  duration  of  this  contract 
specifically  for  these  AE  experiments. 
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3.1.5  Stress-Corrosion  Event  Monitor 

A  conventional  AE  transducer,  in  conjunction  with  a  Dunegan/Endevco  920  Distri¬ 
bution  Analyzer  and  a  Tektronix  604  Monitor,  was  used  to  monitor  the  rate  of  arrival 
and  the  intensity  of  the  stress-corrosion  events.  This  information  allowed  us  to 
categorize  easily  the  state  of  the  C-block  stress  corrosion  AE  source,  as  well  as  provid¬ 
ing  a  useful  reference  standard  against  which  to  measure  the  new  analysis  methods 
under  evaluation. 

3.2  Software 

A  variety  of  software  modules  were  employed  or  developed  to  serve  the  needs  of 
the  experimental  effort.  Software  modules  were  developed  or  used  for  data  acquisi¬ 
tion,  data  validation,  algorithm  identification,  homomorphic  data  analysts,  feature  ex¬ 
traction,  and  pattern  recognition.  Much  of  the  data  analysis  was  carried  out  using  an 
interactive,  interpretive,  high-level  digital  signal  processing  language  available  from 
Signal  Technology  Inc.  (Santa  Barbara,  Calif.)  called  ILS.  “Recipes,”  consisting  of 
DEC  OCL  command  procedures  interpretable  by  ILS,  are  provided  for  all  the  analysis 
carried  out  using  ILS.  In  some  cases,  where  ILS  did  not  provide  appropriate  functional 
modules,  new  ILS  procedures  were  developed  to  provide  the  needed  functionality.  In 
eith<^r  case,  the  modules  used  will  be  functionally  described  below;  source  code,  and 
documentation  for  those  developed  under  this  c  rtract,  can  be  found  in  the  appen¬ 
dices. 

3.2.1  Data  Acquisitioo 

Software  was  developed  to  interface  the  Biomation  hardware  with  the  VAX  11/750 
computer  and  to  provide  some  measure  of  validation  for  the  experimental  data.  Addi¬ 
tional  software  was  developed  to  prepare  the  input  data  in  a  form  compatible  with  the 
high-level  signal  processing  language  (ILS)  which  contains  a  variety  of  software 
modules  for  data  analysis. 

3. 2. 1.1  VAX-Biomation  Interface.  A  data  acquisition/control  software  system  was  put 
together  to  acquire  data  with,  and  control  up  to,  four  Biomation  8100  transient  record¬ 
ers  with  a  VAX  11/750.  The  Biomation  interface  software  is  MACRO-based  and  con¬ 
trols  a  DR-1 1C  parallel  interface  board.  The  interface  permits  VAX-based  software 
control  of  the  Biomation  recorders.  All  of  the  Biomation  front  panel  switches  and  con¬ 
trols  can  be  set  via  software  commands  originating  in  the  VAX.  In  use,  the  Biomation 
transient  recorders  are  operated  in  parallel,  with  one  of  the  recorders  providing  trigger¬ 
ing  synchronization  signals  for  the  others.  The  trigger  delays  were  set  so  that  the  cap¬ 
tured  data  records  would  overlap  slightly.  The  captured  data  records  were  processed  by 
a  VAX  FORTRAN  program  that  carried  out  correlation  calculations  on  the  overlapped 
regions  to  verify  that  the  record  segments  were  correctly  synchronized,  making  ap¬ 
propriate  corrections  if  required.  In  addition,  the  overlapped  regions  were  used  to  ad¬ 
just  the  gain  and  offset  on  each  record  to  compensate  for  record-to-record  variations. 
Finally,  graphics  software  permitted  the  experimenter  to  view  the  captured  records,  or 
any  subset  of  them,  prior  to  proceeding.  Complete  documentation  on  the  system,  in¬ 
cluding  functional  descriptions,  schematics,  and  software  listings  can  be  found  in 
Appendix  A. 


3.2.2  Analysis 

The  basic  approach  in  developing  analysis  software  was  to  carry  out  as  much  of  the 
analysis  as  possible  directly  in  the  ILS  (see  Appendix  B)  digital  signal  processing 
language;  ILS  is  an  expressive  language,  licensed  by  Signal  Technology  Incorporated, 
and  provides  a  flexible  ensemble  of  digital  signal  processing  primitives  built  into  its 
structure.  In  addition,  ILS  provides  many  kinds  of  generalized  support  functions,  in¬ 
cluding  data  editing  capability,  simulation  software,  graphics  software,  and  extensive 
pattern  recognition  software.  ILS  graphics  capability,  for  example,  supports  graphics 
primitives  for  either  interactive  or  hardcopy  display,  along  with  the  associated  support 
functionality. 

Many  complex  functions  can  be  constructed  by  stringing  together  primitive  func¬ 
tional  modules  into  DCL  command  procedure  recipes.  These  recipes  are  executable  ei¬ 
ther  on  a  stand-alone  basis,  or  as  new  primitive  modules  in  more  complex  primitives. 
All  the  recipes  used  in  analyzing  the  data  in  this  report  are  included,  where  appropri¬ 
ate,  with  the  figures  in  Section  4  displaying  processed  data.  Finally,  it  is  relatively  easy 
to  add  entirely  new  modules  to  the  ILS  system. 

Perhaps  the  most  important  feature  of  ILS  is  that  it  is  intended  to  be  used  by 
nonprogrammers;  thus  it  allows  a  researcher  familiar  with  digital  signal  processing 
techniques  to  explore  new  techniques  very  quickly.  This  approach  allowed  considerable 
flexibility,  and  proved  to  be  a  very  effective  way  to  analyze  and  process  quickly  the 
acoustic  emission  data  acquired  under  this  contract. 

3.2.2. 1  Homomorphic  Analysis.  A  software  package  was  developed  under  Phase  I  at 
Rensselaer  Polytechnic  Institute  to  carry  out  simulation  studies  on  the  application  of 
homomorphic  analysis  to  acoustic  emission  signals.  This  package  was  based  on  the 
routines  published  in  the  IEEE  volume  Programs  for  Digital  Signal  Processing.'  After 
some  changes  and  modifications,  the  package  was  installed  in  the  VAX  computer  sys¬ 
tem  and  integrated  with  the  ILS  digital  signal  processing  package  as  the  XCP  module. 
Besides  altering  the  code  to  conform  to  the  ILS  functional  conventions,  two  basic 
changes  were  made  to  the  existing  code.  The  first  change  consisted  of  adding  variable 
exponent  exponential  weighting  capability  to  convert  the  acquired  waveforms  to 
minimum  phase  sequences,  which  are  analytically  more  stable  than  mixed  phase  se¬ 
quences.  (See  Section  2.4  for  a  description  of  the  theory  behind  exponential  weight¬ 
ing.)  The  second  change  involved  adding  bandpass  mapping,  necessary  to  deal  with  the 
signals  produced  by  narrow-band  resonant  transducers.  (See  Section  2.5  for  a  descrip¬ 
tion  of  the  theory  behind  bandpass  mapping.)  The  source  code  and  documentation  for 
the  XCP  ILS  module  can  be  found  in  Appendix  B. 

3.2.2. 2  Feature  Extraction  and  Pattern  Recognition.  A  feature  extraction  and  pattern 
recognition  package  was  put  together  to  investigate  the  separability  of  the  different 
groups  of  waveforms.  The  package  was  based  generally  on  the  pattern  recognition 
software  available  in  ILS.  Since  the  ILS  pattern  recognition  modules  were  originally 
designed  for  speech  and  speaker  recognition,  minor  changes  had  to  be  made  to  some 
of  the  modules  to  make  them  suitable  for  this  application.  The  process  followed  in  us¬ 
ing  the  pattern  recognition  package  to  analyze  acoustic  emission  waveform  data  was 
the  traditional  training/test  set  approach.  A  set  of  waveforms  is  analyzed  and  used  to 
train  the  system;  subsequent  sets  are  tested  against  the  training  set  for  statistical  com¬ 
patibility.  Summarizing  the  process,  each  of  the  waveforms  in  the  group  is  Fourier 
transformed  to  extract  a  number  of  Fourier  components.  These  Fourier  components 
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are  then  processed  to  yield  the  mean  waveform,  the  covariance  matrix,  the  Fischer 
linear  discrimination  matrix,  and  the  eigenvalues  and  eigenvectors  for  each  group  of 
waveforms.  Using  these  intermediate  results,  principal  component  analysis  is  carried 
out  on  sets  of  waveforms.  In  essence,  principal  component  analysis  projects  the  groups 
of  waveforms  onto  a  discriminant  plane  with  maximum  separation  between  the 
groups.  A  scatter  plot  is  then  made  using  the  first  two  principal  components,  and  the 
statistical  properties  of  each  of  the  groups  can  be  computed.  Finally,  the  mean  separa¬ 
tion  between  groups  can  be  expressed  in  statistical  terms  and  used  to  characterize  the 
probability  that  a  given  test  value  is  in  fact  a  member  of  the  training  set.  A  complete 
description  of  the  software  modifications  to  the  ILS  modules  can  be  found  in  Appen- 


Section  4 


EXPERIMENTAL  RESULTS 


The  block  diagram  of  the  experimental  layout,  hardware,  and  associated  electronics 
was  shown  in  Figure  7  of  the  previous  section. 

A  problem  encountered  very  early  in  the  analysis  was  that  our  available  methods 
for  recording  and  for  carrying  out  analog-to-digital  conversion  of  the  data  did  not  pro¬ 
vide  sufficiently  long  record  lengths  for  our  analysis.  Initially,  the  data  acquisition  sys¬ 
tem  utilized  a  single  Biomation  8100  interface  providing  a  pair  of  1024-point-long 
records.  For  sampling  rates  of  10  MHz,  this  system  provided  records  approximately 
0.1  ms  in  duration,  insufficiently  long  for  reverberations  to  decay  in  the  experimental 
fixtures  we  employed.  Accordingly,  several  different  hardware  approaches  for  record¬ 
ing  and  analog-to-digital  conversion  were  explored  to  try  to  identify  a  more  suitable 
means  of  acquiring  the  large  experimental  data  sets  needed  to  confirm  the  results  of 
Phase  I.  The  first  approach  explored  was  to  reduce  the  sampling  frequency  to  2  MHz 
and  to  use  the  Biomation  recorder  in  single-channel  mode,  so  that  a  single  2048-point 
record  is  obtained  instead  of  a  pair  of  1024-point  records.  This  technique  yields  a 
record  duration  of  roughly  1  ms,  which — although  still  too  short  to  record  the  entire 
ringdown  train — provided  enough  data  to  begin  analysis. 

Later  the  record  length  was  increased  fourfold  to  8096  points  by  connecting  four 
Biomation  recorders  in  series,  as  described  above  in  Section  3.1.4.  At  a  sampling  rate 
of  2  MHz,  this  arrangement  gives  a  record  length  of  about  4  ms,  which  was  adequate 
for  subsequent  analysis. 

The  command  procedures  for  generating  the  ILS  plots  shown  in  Figures  9  and  fol¬ 
lowing  are  included  in  the  illustrations.  The  DSP  command  at  the  end  of  most  of  the 
instruction  sets  displays  the  result  in  sampled  data  form. 

4.1  Calibration 

The  breakage  of  Pentel  pencil  leads  was  employed  as  a  calibration  source.  In  the 
initial  phase  of  the  experiments,  the  0.5  mm  HB  leads  were  broken  manually. 

Two  typical  lead-breaking  waveforms  are  shown  in  Figures  9  and  10.  They  are  two 
separate  events  taken  with  two  of  the  wide-band,  NBS-style  conical  transducers  spaced 
about  54  cm  apart  on  the  27  x  5  x  2  in.  aluminum  block.  Pentel  lead  breaks  were 
generated  relatively  close  to  transducer  1  at  the  location  indicated.  Signals  from  trans¬ 
ducer  1  triggered  the  Biomation  to  accept  data.  The  Biomation  was  set  to  the  pretrigger 
mode.  The  data  were  low-passed  at  1  MHz  with  the  Kronhite  filter,  sampled  at  2  MHz, 
for  a  2048-point  record.  The  waveform  shown  in  Figure  9  is  the  output  from  transduc¬ 
er  1,  while  the  waveform  shown  in  Figure  10  is  from  transducer  2. 

The  time  difference  between  the  two  pulses  in  Figures  9  and  10  is  about  200  fis, 
which  is  reasonably  close  to  the  value  of  176  fis  calculated  from  the  geometry  of  the 
setup,  assuming  the  shear  wave  speed  in  aluminum  to  be  3.1  km/s.  The  waveforms 
exhibit  low-frequency  ringdown  due  to  the  reflections  in  the  aluminum  block. 


Homomorphic  deconvolution  was  performed  on  the  waveform  in  Figure  10.  The 
cepstrum  of  the  waveform  is  shown  in  Figure  11.  The  result  of  passing  the  cepstrum 
through  the  low-time  window  w{t): 


1 


wU)  = 


—9.5  fis  <  t  <  9.5  /iS 


(5) 


0  otherwise 


and  then  inverse  transforming  is  shown  in  Figure  12.  An  impulse  type  of  function  is 
recovered,  but  the  result  is  not  very  clean:  the  ringing  in  the  background  is  quite  seri¬ 
ous. 


A  large  improvement  was  brought  about  by  employing  the  technique  of  exponen¬ 
tial  weighting  to  minimize  the  effects  of  reverberation.  The  result  of  applying  exponen¬ 
tial  weighting  to  the  waveform  in  Equation  10  with  a  =  0.9955  and  using  the  low-time 
window  given  by  Equation  5  is  shown  in  Figure  13.  The  recovered  signal  is  much 
sharper  and  cleaner,  and  in  fact  compares  favorably  with  the  theoretical  Pekaris  solu¬ 
tion,  which  is  shown  in  Figure  14. 
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Figure  11.  The  cepstrum  of  the  waveform  in  Figure  10 
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Figure  12.  The  result  in  the  time  domain  of  low>time  pass  filtering  the  cepstrum  of 
the  waveform  in  Figure  10 


4.2  Effect  of  Record  Length 

Another  problem  encountered  was  that  the  ringdown  response  of  the  medium 
I  could  last  substantially  longer  than  the  maximum  record  length,  for  as  long,  in  fact,  as 

several  seconds  under  some  circumstances.  The  record  length  of  a  waveform  deter¬ 
mines  the  fraction  of  ringdown  captured,  and  hence  the  amount  of  leakage  error  in 
the  Fourier  transformation  operation  in  homomorphic  analysis.  One  would  certainly 
expect  the  results  to  improve  substantially  when  the  record  length  is  increased.  As  we 
shall  see,  however,  if  the  effect  of  the  record  length  is  not  significant  enough,  the  im¬ 
provement  brought  about  by  longer  record  length  could  easily  be  masked  by  other 
sources  of  errors. 


The  effect  of  the  record  length  is  demonstrated  in  the  results  shown  in  the  se¬ 
quence  of  Figures  15-26.  Figures  15-18  are  four  lead-breaking  waveforms  captured 
with  a  conical  transducer,  the  source  being  located  at  distances  45,  45,  30,  and  15  cm, 
respectively,  from  the  transducer.  The  signals  were  low-pass  filtered  at  1  MHz,  sam¬ 
pled  at  2  MHz,  and  are  8096  points  in  length.  Each  of  these  waveforms  was 
transformed  to  the  cepstral  domain,  then  low-time  filtered  with  the  wit): 


—4.5  fis  <  t  <  4.5  fis 


0  otherwise 


(6) 
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Figure  13.  The  improved  result  obtained  through  the  use  of  exponential  weighting; 
a  -  0.9955 


and  finally  inverse-transformed  to  the  time  domain.  The  value  of  a  used  in  exponen¬ 
tial  weighting  is  0.9985.  Figures  19  through  22  show  the  results  obtained  in  this  way 
using  all  of  the  8096  samples  of  the  input  waveforms  in  Figures  15  through  18.  In 
comparison,  the  results  obtained  using  only  the  first  2048  samples  of  the  input 
waveforms  are  shown  in  Figures  23  through  26.  Identical  values  of  a  and  wit)  are 
used  in  both  cases.  In  comparing  the  two  sets  of  results,  although  Figure  19  is  much 
cleaner  than  Figure  23,  the  improvement  caused  by  the  longer  record  length  is  less 
obvious  in  the  other  three  cases.  From  these  results  it  may  be  concluded  that  the  leak¬ 
age  error  in  Fourier  transformation  in  going  from  4.096  ms  (corresponding  to  8192 
points  at  2  MHz  sampling  rate)  to  1.024  ms  is  not  significant  compared  to  the  other 
sources  of  errors  in  the  experiments  and  analysis.  The  lack  of  improvement  with 
longer  record  length  may  also  be  due  partially  to  some  transient  instabilities  in  the  gain 
and  noise  characteristics  of  the  four  Biomation  recorders.  At  any  rate,  the  record 
length  of  1.024  ms  provided  by  one  Biomation  recorder  at  a  sampling  rate  of  2  MHz 
appeared  to  be  long  enough  for  our  purpose.  The  waveforms  used  in  some  of  our  later 
analyses  were  1.024  ms  in  length. 
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Figure  14.  The  theoretical  seismic  surface  pulse  according  to  Pekaris 


4.3  Dependence  on  the  Width  of  the  Cepstral  Filter 

The  ideal  width  of  the  cepstral  filter  used  to  separate  out  the  signal  impulses  from 
the  recorded  waveforms  was  determined  empirically.  The  results  indicated  that  in  most 
cases  a  filter  centered  around  the  origin  in  the  cepstral  domain  with  a  width  of  about 
4.5  fjiS  worked  best.  Two  typical  sets  of  results  are  illustrated  in  Figures  27-32.  Figures 
27-29  show  the  results  of  low-passing  the  cepstrum  of  the  recorded  waveform  in  Fig¬ 
ure  10  using  a  =  0.999  in  exponential  weighting  and  cepstral  filters  of  widths  ±2.5, 
±4.5,  and  ±9.5  fis,  respectively.  Figures  30-32  are  the  results  of  applying  the  same 
procedure  to  the  waveform  in  Figure  15  using  a  =  0.9985  in  exponential  weighting.  In 
both  cases,  the  cepstral  filter  with  a  width  of  ±4.5  /xs  yielded  the  best  results. 

4.4  Dependence  on  Alpha 

In  the  homomorphic  analysis  of  single-event  waveforms,  exponential  weighting 
converts  the  original  time  sequences  hp{n)  and  hrin)  to  the  new  ones  a” hpin)  and 
a"hr(n).  The  value  of  a  determines  the  distributions  of  the  cepstra  of  the  two  ex¬ 
ponentially  weighted  time  sequences,  and  hence  the  results  of  the  homomorphic 
deconvolution  depend  on  the  value  of  a  used  in  the  exponential  weighting.  To  study 
the  extent  of  this  dependence,  homomorphic  analysis  was  performed  on  the  waveform 
in  Figure  15  with  the  window  wf/)  in  Equation  6  but  with  different  values  of  a  in  the 
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Figure  18. 
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A  lead'break  waveform  recorded  by  an  NBS  conical  transducer  located 
30  cm  from  the  source 
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A  lead-break  waveform  recorded  by  an  NBS  conical  transducer  located 
15  cm  from  the  source 
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Figure  19.  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  15,  using  ail  of  the  8096  samples  as  input;  a  —  0.9985 


exponential  weighting.  The  results  corresponding  to  o  -  0.9985,  0.998,  and  0.9975  are 
shown  in  Figures  33,  34,  and  35,  respectively.  The  results  indicate  that  even  a  slight 
change  in  the  value  of  a  affects  the  results  significantly. 

4.5  Fourier  Deconvolution 

Initial  efforts  to  use  the  transfer  function  obtained  in  the  homomorphic  deconvolu¬ 
tion  of  one  Pentel  lead-breaking  waveform  to  deconvolve  other  Pentel  lead-breaking 
waveforms  were  met  with  little  success.  To  understand  the  difficulties,  the  repeatability 
of  the  manually  performed  lead-breaking  waveforms  was  examined  in  detail.  It  was 
found  that  the  lead-breaking  signals  depend  critically  on  (1)  the  precise  angle  the  lead 
makes  with  the  surface  and  (2)  possible  secondary  Pentel  tip  impact  after  the  lead 
break. 

These  findings  were  illustrated  in  Figures  36-38,  which  are  manually  performed 
lead-breaking  waveforms  captured  with  a  conical  transducer.  Figures  36  and  37  illus¬ 
trate  the  dependence  of  the  signal  on  the  angle  between  the  lead  and  the  surface.  The 
waveform  in  Figure  36  was  generated  with  an  angle  of  about  65°,  and  the  one  in  Fig¬ 
ure  37  with  an  angle  of  about  40°.  The  waveforms  are  markedly  different.  It  appears 
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Figure  22.  The  result  of  low>time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  18,  using  all  of  the  8096  samples  as  input;  a  ==  0.9985 


that  the  angle  determines  the  amount  of  shear  and  transverse  waves  generated  in  each 
event.  A  shallow-angle  (<45°)  event  contains  relatively  more  low-frequency  com¬ 
ponents  and  is  relatively  featureless.  On  the  other  hand,  a  steep  angle  (~60°-70°) 
produces  strongly  featured  signals  with  much  modulation.  A  Pentel  lead  break  with 
secondary  Pentel  tip  impact  at  large  angle  is  shown  in  Figure  38.  The  secondary  impact 
generated  additional  AE  activity  in  the  midregion  of  the  waveform. 

In  an  effort  to  produce  repeatable  lead-breaking  signals,  a  fixture  for  breaking  leads 
in  a  well-controlled  manner  was  acquired.  The  fixture  allowed  precise  control  over  the 
angle  of  the  lead  during  breakage  and  effectively  eliminated  problems  with  periodic 
secondary  tip  contact.  Further  studies  on  the  feasibility  of  Fourier  deconvolution  were 
carried  out  using  waveforms  generated  with  the  help  of  the  fixture.  It  was  found  that 
the  waveforms  were  highly  repeatable.  The  repeatability  is  demonstrated  in  Figures  39 
and  40,  where  200  points  of  two  such  waveforms  are  shown  near  the  initial  impulses. 
These  waveforms  were  recorded  by  a  conical  transducer  and  sampled  at  2  MHz.  Fig¬ 
ures  41  and  42  show  another  200  points  in  the  later  part  of  the  two  waveforms.  Minor 
differences  between  the  two  waveforms  are  visible  after  a  considerable  delay,  but  the 
overall  similarity  is  still  very  impressive. 

Fourier  deconvolution  using  these  mechanically  generated  waveforms  yielded  en¬ 
couraging  results.  Figure  43  shows  the  impulse  response  obtained  by  the  following 
procedure; 

•  Averaging  33  of  such  waveforms  (Pentel  lead  breaking  by  a  mechanical  device 
recorded  by  a  conical  transducer),  each  of  the  waveforms  being  1024  points  in 
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Figure  23.  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  IS,  using  only  the  first  2048  samples  as  input;  a  =  0.9985 


•  Calculating  the  cepstrum  of  the  averaged  waveform,  a  =  0.999 

•  High-pass  filtering  the  cepstrum,  i.e.,  setting  the  low-time  portion  of  the  cep¬ 
strum  to  zero  for  |r|  <  4.5  fis 

•  Inverse-transforming  the  cepstrum  to  the  time  domain 

The  result  of  deconvolving  one  of  the  33  waveforms  used  in  averaging,  the  one 
shown  in  Figure  44,  by  the  impulse  response  in  Figure  43  is  illustrated  in  Figure  45. 
The  result  was  improved  by  removing  the  high-frequency  noise  with  the  elliptical  fre¬ 
quency  filter  shown  in  Figure  46  with  —60  dB  cut-olf  at  700  kHz;  the  filtered  result  is 
show  in  Figure  47. 

The  deconvolution  results  shown  in  Figures  45  and  47  were  obtained  under  close 
to  ideal  conditions: 

•  Obtaining  highly  repeatable  waveforms  by  breaking  Pentel  leads  with  a  mechani¬ 
cal  device 
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Figure  26.  The  result  of  low>time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  18,  using  only  the  first  2048  samples  as  input;  a  =  0.9985 


•  Recording  the  waveforms  with  a  wide-band  conical  transducer 

•  Using  a  relatively  short  record  of  1024  points,  where  the  similarity  between 
waveforms  is  extremely  good 

In  Figures  45  and  47,  the  presence  of  the  impulse  is  obvious,  yet  the  quality  of  the 
pulse  shape  does  not  appear  to  be  good  enough  for  identification  purposes,  even  after 
low-pass  filtering  to  remove  high-frequency  noise.  These  results  indicate  that  Fourier 
deconvolution  is  too  sensitive  to  the  minor  differences  in  the  recorded  waveforms  to 
be  of  any  value  for  event  identification  in  practical  situations. 

4.6  Effects  of  a  Limited-Frequency  Band 

Realistic  acoustic  emission  events  were  generated  using  aluminum  C-blocks  subject 
to  stress  and  a  salt  solution.  It  was  found  that  the  amplitude  of  stress-corrosion  events 
generated  in  this  way  was  about  60  dB  lower  than  the  amplitude  from  Pentel  lead 
breaks;  this  finding  can  be  seen  clearly  in  the  amplitude  plot  shown  in  Figure  48.  Un¬ 
fortunately,  since  the  sensitivities  of  the  conical  transducers  are  much  lower  than 
those  of  the  resonant  transducers,  the  stress-corrosion  signals  could  not  be  detected  by 
the  conical  transducers  at  all.  For  this  reason,  resonance  transducers  were  used  for 
detection. 

The  data  taken  with  the  resonance  transducers  are  band-limited  between  100  kHz 
and  300  kHz.  The  power  spectrum  of  one  of  the  stress-corrosion  waveforms  captured 
with  a  resonance  transducer  is  shown  in  Figure  49.  As  mentioned  in  Section  2.5,  such 
data  require  bandpass  mapping  to  suppress  the  effects  of  noise  propagation  in  the 
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Figure  28.  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  10,  a  =  0.999,  width  of  cepstral  filter  =  ±4.5  /lis 


Figure  51  indicates  that  the  band-limited  nature  of  the  waveforms  recorded  using 
resonance  transducers  gives  rise  to  severe  distortions  in  the  impulses  obtained  by  low- 
pass  filtering  the  cepstrum  of  the  averaged  waveform  in  Figure  50.  Bandpass  mapping 
cleans  up  the  result  considerably  by  suppressing  the  error  propagation  in  the  unoccu¬ 
pied  frequency  regions,  but  severe  ringing  and  broadening  is  still  present  as  shown  in 
Figure  52.  This  ringing  and  broadening  is  due  to  intrinsic  limitations  of  the  band- 
limited  nature  Oi  the  data — resulting  from  the  use  of  resonance  transducers — and 
therefore  cannot  be  removed  by  signal  processing. 

4.7  Pattern  Recognition 

The  quality  of  the  homomorphically  deconvolved  signals  obtained  above  does  not 
give  one  much  hope  for  accomplishing  source  identification.  This  observation  was 
confirmed  in  the  following  pattern  recognition  analysis  using  two  sets  of  waveforms 
from  different  sources,  with  nine  waveforms  in  each  set.  In  one  set  are  waveforms 
generated  by  rubbing  sandpaper  on  the  edge  of  one  end  of  the  aluminum  block,  while 
the  ones  in  the  other  set  were  generated  by  rubbing  the  sandpaper  on  the  face  of  the 
same  end  of  the  block.  The  waveforms  were  sampled  at  2  MHz  and  are  2048  points 
in  length.  Each  of  the  waveforms  was  homomorphically  deconvolved  to  recover  the 
low-time  component  in  the  cepstral  domain,  using  a  weighting  of  0.9985  and  a  cepstral 
cutoff  time  of  ±4.5  fis.  These  deconvolved  signals  were  Fourier  analyzed,  and  the  20 
Fourier  coefficients  from  976.56  Hz  to  594.7  kHz  at  31.25  kHz  intervals  were  extract¬ 
ed.  Principal  component  analysis*’  was  performed  on  these  20  Fourier  coefficients  of 
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Figure  29.  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  10,  a  =  0.999,  width  of  cepstral  filter  =  ±9.5  ixs 


the  two  sets  of  9  waveforms  in  each  set  using  the  software  in  the  ILS  signal  analysis 
package.  The  first  two  principal  components  are  plotted  in  Figure  53.  (See  Appendix  C 
for  a  description  on  pattern  recognition  and  examples  of  command  procedures  for  Fig¬ 
ures  53  and  following.)  The  results  in  Figure  53  indicate  the  two  sets  of  waveforms 
cannot  be  distinguished  from  each  other.  In  comparison,  the  first  two  principal  com¬ 
ponents  of  the  raw  data  without  cepstral  filtering  are  plotted  in  Figure  54.  The  two  sets 
of  waveforms  become  less  separated  after  the  cepstral  filtering. 

Similar  results  were  obtained  when  the  above  procedure  was  repeated  with  an  addi¬ 
tional  set  of  nine  lead-breaking  waveforms.  The  plots  of  the  first  two  principal  com¬ 
ponents  with  and  without  cepstral  filtering  are  illustrated  in  Figures  55  and  56.  Again, 
cepstral  filtering  failed  to  improve  the  separability  of  the  three  sets  of  waveforms. 
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Figure  30.  The  result  of  low*time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  15,  a  =  0.9985,  width  of  cepstral  filter  =  ±2.5  ^is 
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Figure  31.  The  result  of  low>time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  15,  a  ^  0.9985,  width  of  cepstral  filter  —  ±4.5  us 
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Figure  32.  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  15,  a  =  0,9985,  width  of  cepstral  filter  =  ±9.5  fj.s 
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Figure  33.  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  waveform 
in  Figure  15,  a  —  0.9985 


cm  2i4a 

flL  12$$ 
OSf 


'■ID  •  .••21)46  -ec 


(HKASURtO  WAVCrOM  WO  PIOS  NO.) 


Figure  37.  A  typical  small-angle  Pentei  lead-break  waveform;  the  angle  between 
the  lead  and  the  surface  is  40° 


Figure  38.  A  Pentel  lead-break  waveform  with  secondary  Pentel  tip  impact  at  a 
large  angle 
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Figure  39.  A  200-point  sample  of  a  Bxture-generated  lead-break  waveform  near  the 
initial  impulse 


Figure  41.  Another  200-point  sampie  of 
part  of  the  waveform 
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Figure  43.  The  impulse  response  obtained  by  high-time  pass  filtering  the  cepstrum 
of  the  average  of  33  mechanically  generated  lead-break  waveforms;  a  = 
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The  result  of  deconvolving  the  waveform  in  Figure  44  by  the  impulse 
response  in  Figure  43 
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The  resuit  of  iiitering  the  waveform  in  Figure  45  with  the  elliptical  Alter 
in  Figure  46 
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Relative  amplitude  of  the  stress  corrosion  events  and  the  lead-breaking 
events 
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Figure  49.  Power  spectrum  of  a  stress-corrosion  event  captured  by  a  resonance 
transducer 
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Figure  50.  The  average  of  31  stress-corrosion  waveforms  with  their  initial  impulses 
lined  up  at  the  same  location;  the  waveforms  were  captured  with  a  reso¬ 
nance  transducer 
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Figure  52.  The  result  of  low-time  pass  filtering  the  cepstrum  of  the  averaged 
waveform  in  Figure  50;  band-pass  mapping  between  100  kHz  and 
300  kHz 
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Figure  S3.  Scatter  plots  of  the  first  two  principal  components  of  the  frequency  com¬ 
ponents  of  two  sets  of  events  after  low-time  pass  filtering  in  the  cepstral 
domain.  The  two  sets  of  events  were  generated  by  (a)  rubbing  a  sandpa¬ 
per  on  the  edge  of  one  end  of  the  aluminum  block  and  (b)  nibbing  the 
sandpaper  on  the  facet  of  the  same  end  of  the  block 
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Figure  54.  Scatter  plots  of  the  first  two  principal  components  of  the  frequency  com* 
ponents  of  the  two  sets  of  events  in  Figure  53  without  filtering  in  the 
cepstral  domain 


Figure  55.  Scatter  plots  of  the  first  two  principal  components  of  the  frequency  com¬ 
ponents  of  three  sets  of  events  after  low-time  pass  filtering  in  the  cep- 
stral  domain;  the  three  sets  of  events  were  generated  by  (a)  rubbing  a 
sandpaper  on  the  edge  of  one  end  of  the  aluminum  block,  (b)  rubbing 
the  sandpaper  on  the  facet  of  the  same  end  of  the  block,  and  (c)  lead¬ 
breaking 


Section  5 


CONCLUSION 


The  key  results  obtained  in  this  study  are  summarized  below: 

•  Adaptive  homomorphic  deconvolution  seems  to  work  relatively  well  for  AE  sig¬ 
nals  captured  with  a  wide-band  flat  response  transducer  with  good  S/N.  This 
technique  is  able  to  clean  up  multipath  and  multimode  reverberation  and  ring- 
down  effects  rather  well.  The  technique  does,  however,  have  some  unpleasant  ad 
hoc  elements,  in  that  several  key  parameters  cannot  be  determined  a  priori. 
Specifically,  the  decay  exponent,  needed  for  exponential  weighting  to  make  the 
sequence  a  well-behaved  minimum  phase  sequence,  is  a  problem,  as  are  the 
linear  filter  coefficients  needed  for  filtering  the  cepstral  waveform.  The  appropri¬ 
ate  constants  can  quickly  be  found  with  a  little  experimentation,  but  the  process 
is  non-unique. 

•  Adaptive  homomorphic  deconvolution  seems  to  work  to  some  degree  with  a 
conventional  resonant  transducer,  particularly  when  bandpass  mapping  is  em¬ 
ployed.  The  success  of  the  deconvolution  is  severely  limited,  however,  because 
of  the  limited  signal  bandwidth  information  provided  by  the  resonant  transducer. 
As  can  be  seen  in  plots  shown  in  Section  4.6,  although  pulse  compression  has 
been  achieved,  the  reconstructed  signal  suffers  from  serious  ringing  effects  due 
to  the  very  limited  transducer  bandwidth. 

•  Adaptive  homomorphic  deconvolution,  implemented  by  high-pass  filtering  the 
averaged  cepstral  domain  signal,  seems  to  be  effective  in  providing  approximate 
estimates  of  the  workpiece  transfer  function,  particularly  when  the  original  signal 
is  obtained  with  a  wide-band  transducer.  As  noted  above,  the  technique  is  not 
effective  when  narrow-band  resonant  transducers  are  employed,  because  of  the 
limited  illumination  of  the  cepstral  domain  afforded  by  those  transducers.  Un¬ 
fortunately,  the  wide-band  transducers  employed  in  this  study  are  not  suitable 
for  use  in  practical  applications  because  of  their  fragility  and  severely  limited  sen¬ 
sitivity.  Their  sensitivity  is  not  high  enough  to  be  used  to  obtain  signals  from 
typical,  real  AE  sources;  specifically,  they  are  too  insensitive  for  use  with  alumi¬ 
num  stress  corrosion  V-block  sources,  or  most  practical  AE  sources  in  general. 

•  The  final  step  in  employing  adaptive  homomorphic  deconvolution  to  analyze 
acoustic  emission  signals  is  to  employ  ordinary  Fourier  deconvolution  to  elimi¬ 
nate  the  effects  of  transducer  ringdown,  multipath,  etc.  Fourier  deconvolution, 
when  used  to  eliminate  the  workpiece  transfer  function  effects,  is  very  sensitive 
to  the  precise  details  of  the  transfer  function.  In  practical  terms,  this  means  that 
the  transfer  functions  must  be  precisely  repeatable  from  pulse  to  pulse.  In  fact, 
the  data  obtained  in  this  study  leads  one  to  believe  that  there  are  substantial 
variations  in  the  transfer  function  obtained  from  both  artificial  and  natural  AE 
sources,  particularly  when  relatively  long  records  of  the  pulse  ringdown  signals 
are  obtained.  The  transfer  functions  are  more  repeatable  when  relatively  short 
records  are  used.  Preliminary  analysis  suggests  that  slight  variations  in  the  direc¬ 
tion  of  the  signal  propagation  may  be  responsible  for  these  effects.  A  slight  varia¬ 
tion  in  the  direction  of  propagation,  or  in  the  mix  of  longitudinal  and  shear 
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waves,  will  result  in  changes  in  the  workpiece  ringdown  response  function;  these 
changes  get  larger  and  larger  as  the  ringdown  proceeds.  Thus,  the  later  segments 
of  long  records  will  suffer  progressively  increasing  distortion  from  inaccurate 
deconvolution.  Natural  AE  signals,  since  they  emanate  from  evolving  cracks 
where  the  crack  propagation  follows  natural  grain  boundaries,  may  be  expected 
to  offer  some  variation  in  the  strength  and  direction  of  the  emitted  shear  and 
longitudinal  waves.  Surprisingly,  even  AE  signals  carefully  generated  via  artificial 
means  (Pental  lead  breaks)  seem  to  vary  in  the  precise  details  of  their  ringdown 
responses  over  relatively  long  record  lengths.  The  net  result  is  that  only  rela¬ 
tively  short  (~2048  points)  record  lengths  can  be  used  successfully,  either  with 
natural  or  artificial  AE  signal  sources. 

•  Pattern  recognition,  employed  as  a  means  of  identifying  and  characterizing  AE 
sources  through  differences  in  the  microstructure  of  their  waveforms,  does  not 
seem  to  be  particularly  effective,  either  when  applied  to  the  original  waveforms, 
or  when  applied  to  waveforms  restored  using  adaptive  homomorphic  deconvolu¬ 
tion.  The  original  waveforms  are  heavily  contaminated  with  the  workpiece 
transfer  function  effects,  and  microstructure  effects  are  difficult  to  separate.  The 
reconstructed  waveforms  suffered  from  ringing  distortion  because  of  the  limited 
bandwidth  of  the  transducers  employed.  They  were  also  very  difficult  to  separate. 

In  conclusion,  the  analytical  techniques  explored  in  this  research  provide  several 
new  options  in  processing  acoustic  emission  signals.  Specifically,  they  provide  the 
means  to  obtain  and  employ  workpiece  transfer  functions  in  deriving  the  underlying 
acoustic  emission  waveforms.  Successful  application  of  these  techniques  requires  the 
use  of  wide-band  transducers  that  are  faithful  to  the  original  waveforms  and  can  only 
be  employed  on  relatively  short  records.  Any  practical  application  must  await  the 
development  of  wide-band  accurate  transducers  with  much  improved  sensitivity  over 
currently  available  transducers. 
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Appendix  B 

ILS  SOFTWARE  MODULES 


The  ILS  Interactive  Laboratory  System^"  interactive  digital  signal  processing 
software  package  (available  through  software  leasing  agreements  from  Signal  Technol¬ 
ogy  Inc.,  Santa  Barbara,  California)  was  used  as  the  primary  software  tool  in  develop¬ 
ing  the  signal  processing  software  used  in  this  study.  The  basic  ILS  package  supplies  a 
very  flexible  and  powerful  set  of  modular  software  components,  implemented  in  an  in¬ 
teractive  environment  featuring  considerable  graphics  feedback.  The  ILS  software  pro¬ 
vides  many  generic  digital  signal  processing  modules  useful  in  a  wide  variety  of  appli¬ 
cation  areas.  When  we  found  that  ILS  did  not  provirle  the  necessary  functionality 
needed  in  the  acoustic  emission  application,  new  ILS  modules  were  written,  or 
modifications  were  made  to  existing  ILS  modules  to  implement  the  new  functionality. 
These  software  modules  and  documentation  art  included  on  the  following  pages  of 
this  Appendix. 
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Existing  ILS  commands  as  supplied  by  STI; 


ADF  ARITHMETIC  FUNCTIONS  ON  INTEGER  DATA 
AFP  3-D  AREA  FUNCTION  PLOT  OP  ANALYSIS  DATA 
ANL  LINEAR  PREDICTION  ANALYSIS  USING  AUTOCORRELATION, 
COVARIANCE,  OR  BURG'S  METHOD 
API  INVERSE  FILTER  ANALYSIS  USING  AUTOCORRELATION 

METHOD  WITH  PITCH  DETECTION  USING  CEPSTRAL  METHOD 
ASG  ASSIGN  LOGICAL  UNIT  NUMBER  FOR  KEYBOARD 
AVG  RECORD  AVERAGING  (MOVING  OR  EXPONENTIAL) 

INPUT  IS  PRIMARY  PILE,  OUTPUT  SECONDARY 
BOP  BINARY  ARITHMETIC  OPERATIONS  ON  FLOATING  POINT 
DATA;  OR  COMBINE  TWO  REAL  VECTORS  INTO  ONE 
COMPLEX  VECTOR  -  TWO  PILES  INPUT,  ONE  OUTPUT 
BPA  PATTERN  MATCHING  OF  TEST  AND  REFERENCE  DATA 

SUPPORTS  DYNAMIC  PROGRAMMING  FOR  APPLICATIONS 
SUCH  AS  WORD  RECOGNITION 

CEP  CEPSTRUM  DISPLAY  OF  VARIABLE  LENGTH  FRAMES 
CLA  CURSOR  LABELING  OF  PERIODIC  SAMPLED  DATA 
CNV  CONVOLUTION  OP  A  LONG  DURATION  TIME  SERIES 
WITH  A  FILTER  IMPULSE  RESPONSE  USING  SHORT 
OVERLAPPED  FPT'S 

COR  CORRELATION  OP  TWO  LONG  DURATION  TIME  SERIES 
USING  SHORT  OVERLAPPED  FFT'S 

CST  COMPUTATION  OP  SAMPLED  DATA  STATISTICS;  PRINTS 

LOCATION  OF  PEAK  VALUE 

CTX  EXAMINE  OR  CHANGE  CONTEXT,  WHICH  IS  THE  NUMBER 

OF  SAMPLED  DATA  POINTS  PER  FRAME  OF  DATA 
CUR  X-Y  CURSOR,  USED  IN  CONJUNCTION  WITH  DISPLAY 
DPM  DISPLAY  OF  ANALYSIS  DATA;  SUPPORTS  VARIABLE 
FRAME  SIZES 

DRE  DISPLAY  REAL  OR  COMPLEX  VECTORS  FROM  RECORD  FILES 

DSP  DISPLAY  OF  TIME  SERIES;  GRIDS  AND  FRAME  LOCATIONS 

EFI  ELLIPTIC,  BUTTERWORTH  AND  CHEBYCHEV  FILTER  DESIGN 

FDI  FREQUENCY  SPECTRUM  DISPLAY,  PROM  FFT  ON  SAMPLED 

DATA  (MAXIMUM  OF  512  POINTS) 

FFT  FAST  FOURIER  TRANSFORM  OF  A  TIME  SERIES;  SUPPORTS 
CIRCULAR  SHIFTS  OP  THE  DATA 
FIL  SPECIFY,  CREATE,  DELETE  BINARY  DATA  FILES 

FLT  CONVOLUTION  OP  A  TIME  SERIES  WITH  DIGITAL 

FILTERS;  SUPPORTS  INTEGER  OR  FLOATING  POINT 
DATA  FILES 

FPL  FREQUENCY  PLOT  OF  SMOOTH  SPECTRA  FROM 
ANALYSIS  FILE 

FTR  FORMANT  TRACKING  OF  SPEECH  DATA;  STORES 

BANDWIDTH  INFORMATION  AND  SUPPORTS  VARIABLE 
FRAME  SIZES 

GRD  PLOT  A  GRID  FOR  FREQUENCY  OR  AN  AXIS  FOR  TIME 
HHH  HELP  WITH  ILS  COMMANDS 

HIS  READ  DATA  FROM  RECORD  FILES  AND  PLOT  A 
HISTOGRAM  ON  THE  SCREEN 
IDC  ENTER  DATA  INTO  COMMON  BLOCK 


IFL  IDEAL  FILTER  DESIGN;  SUPPORTS  CIRCULAR  SHIFT 
OF  IMPULSE  RESPONSE 

ILS  CREATE  AND  INITIALIZE  USER  COMMON  FILE 
INA  INITIALIZATION  OF  SAMPLED  DATA  FILE  HEADER; 

SUPPORTS  SPECIFICATION  OF  LOGICAL  END  OF  FILE 
WHICH  IS  USED  IF  STARTING  FRAME  AND  NUMBER 
OF  FRAMES»1  IN  SUBSEQUENT  PROGRAMS 
LBA  LABEL  A  SAMPLED  DATA  SEGMENT 
LBF  POINT  TO  A  LABEL  FILE 
LCM  LIST  ILS  COMMON  (ALL  OR  PART) 

LFI  DESIGNS  LINEAR  PHASE  FINITE  IMPULSE  RESPONSE 

(FIR)  FILTERS  USING  THE  REMEX  EXCHANGE  ALGORITHM 
LLA  LISTS  LABEL  FILE 

LRE  LISTING  OF  FLOATING  POINT  DATA;  INDEXES  TIME 

SERIES  IN  SECONDS  AND  SPECTRAL  DATA  IN  HERTZ 
LSN  DIGITAL-TO-ANALOG  CONVERSION  *NOT  IMPLEMENTED* 

MDF  MODIFY  VALUES  OF  DATA  POINTS  IN  FILES 

WILL  MODIFY  FILE  HEADERS  AND  RECORD  HEADERS 
MDX  MULTIPLEXING  AND  DEMULTIPLEXING  OF  MULTI¬ 
CHANNEL  SAMPLED  DATA  FILES 

MRE  MANIPULATE  RECORD  FILES  BY  EXTRACTING  RECORDS, 

ITEMS  AND  ELEMENTS 

MVF  DATA  MOVING  WITHIN  A  FILE;  SUPPORTS  ANALYSIS 
FRAMES  AND  ZEROING  OF  ORIGINAL  DATA 
NSI  SIMULATION  OF  NOISY  DATA  BY  ADDING  NOISE  TO  A 

SIGNAL;  GENERATES  PSEUDORANDOM  NOISE 
OPN  ALLOCATES  AND  OPENS  RECORD  FILES 

PAN  PITCH  SYNCHRONOUS  ANALYSIS  OF  SAMPLED  DATA 

PCO  PRINCIPAL  COMPONENTS  ANALYSIS 

PLR  SCATTER-PLOT  OF  DATA  ELEMENTS  IN  RECORD  FILES 
PNS  PITCH  SYNCHRONOUS  SYNTHESIS;  INCLUDES  IMPROVED 
ALGORITHMS  WITH  GLOTTAL  PULSE  OR  GAUSSIAN 
NOISE  EXITATION 

PRT  PRINT  PROM  BINARY  FILES  (OR  HEADER) 

QUR  QUEUEING  OF  FEATURES  IN  TIME  SERIES  BASED  ON 

LABEL  INFORMATION  FOR  PATTERN  RECOGNITION. 
INCLUDES  GENERATION  OF  FEATURE  MATRICES  FOR 
APPLICATIONS  SUCH  AS  PATTERN  RECOGNITION. 

RAN  RESIDUE  ANALYSIS  USING  COEFFICIENTS  PROM  AUTO¬ 
CORRELATION  OR  COVARIANCE  ANALYSIS 
REC  ANALOG-TO-DIGITAL  CONVERSION  INTO  A  FILE 
IS  NOT  IMPLEMENTED 

EQUIVALENT  HARDWARE  AND  SOFTWARE  EXISTS 
FOR  SOME  APPLICATIONS 

RSO  ROOT  SOLVING  FOR  SPECTRAL  RESONANCES;  STORES 
REAL  ROOTS  FOR  HIGH  QUALITY  SPEECH  SYNTHESIS 
RVR  REVERSE  THE  ORDER  OP  SAMPLED  DATA  POINTS 
SDE  SPECTRAL-DENSITY  ESTIMATION  COMMAND 

(CROSS  OR  AUTO) 

SDI  3-DIMENSIONAL  DISPLAY  OF  SPECTRA;  SUPPORTS 

INTEGER  TIME  SERIES,  FLOATING  POINT  TIME 
SERIES  OR  FLOATING  POINT  SPECTRA  INPUT  PILES 
SGM  SPECTROGRAM  DISPLAY  OF  ANALYSIS  DATA;  SUPPORTS 

VARIABLE  FRAME  SIZES  AND  ERASES  BEFORE  DISPLAY 


SIF  FUNDAMENTAL  FREQUENCY  EXTRACTION  USING  SIFT 
ALGORITHM 

SME  CALCULATES  STATISTICS  OF  DATA  IN  FEATURE  RECORDS 
SNS  SYNTHESIZE  SAMPLED  DATA  FROM  ANALYSIS  DATA 
SPL  3-DIMENSIONAL  SPECTRAL  PLOT  OF  ANALYSIS  DATA; 

SUPPORTS  VARIABLE  FRAME  SIZES 
SRE  STORE  RECORDS  INTO  SECONDARY  FILE  FROM 

KEYBOARD  INPUT,  PRIMARY  SAMPLED  DATA  FILE  OR 
PRIMARY  ANALYSIS  FILE 

SSP  FREQUENCY  SPECTRUM  DISPLAY  FROM  FFT  OF  INVERSE 
FILTER  COEFFICIENTS 

TBL  SET  UP  DIRECTORY  TABLE  FOR  ILS  FILE  SYSTEM 
TFU  PROGRAM  FOR  CREATING  TEST  DATA  IN  A  FILE 
(FILE  MUST  BE  CREATED  FIRST) 

TLA  COPY  LABELS  MEETING  A  GIVEN  SPECIFICATION  INTO 
A  SECONDARY  LABEL  FILE 

TRE  TRANSFER  OF  FLOATING  POINT  DATA;  SUPPORTS 
CHANGES  IN  RECORD  SIZE 

TRF  DATA  TRANSFER  BETWEEN  FILES;  SUPPORTS  ANALYSIS 
FRAMES 

TRM  SETTING  OF  TERMINAL  CHARACTERISTICS;  SUPPORTS 
NUMBER  OP  GRAPHIC  INPUT  TERMINATORS  AND  THE 
HP2648  TERMINAL  IN  NATIVE  MODE.  ALSO  PRINTS 
SYSTEM  TIME  AMD  DATE  AMD  WILL  LABEL  A  PLOT. 

WILL  CLEAR  SCREEN  OR  PUT  TERMINAL  IN  ALPHA  MODE 
TSI  TEST  SIGNAL,  GENERATE  SAMPLED  DATA,  LABEL 
OR  RECORD  FILES. 

TTL  TRANSFERS  MARKED  SECTION  OP  SAMPLED  DATA  TO 
THE  SECONDARY  PILE  WITH  OPTIONAL  LABELING 
UOP  UNARY  OPERATIONS,  PHASE  UNWRAPPING  OF  FFT  DATA; 
USES  A  MUCH  IMPROVED  ALGORITHM.  MANY  OTHER 
MANIPULATIONS  OF  FLOATING  POINT  DATA 
VDI  VARIABLE  DISTANCE  THRESHOLD  EVALUATION 

VER  VERIFY  HEADER  BLOCK  IN  SAMPLED  DATA  OR 
ANALYSIS  PILE 

VTR  PLOT  A  VOCAL  TRACT  PROM  SECONDARY  ANALYSIS 
FILE  OR  USER'S  COMMON 

XPA  EXPAND,  INTERPOLATE  OR  DOWNSAMPLE  PRIMARY 

SAMPLED  DATA  FILE  FOR  HIGHER  OR  LOWER  SAMPLING 
FREQUENCIES,  BY  INTERLEAVING  ZEROS  WITH  THE 
ORIGINAL  DATA  OR  SKIPPING  DATA 
XTR  READ  DATA  FROM  PRIMARY  RECORD  FILE  AND  COMPUTE 
MAXIMUM  OR  MINIMUM  VALUES 


******************************************************************* 

New  ILS  commands  developed  and  used  in  the  course  of  this 
research. 

mt****  It**  ***********  ************************  **1111  hit***************** 


$  APP  -  MERGES  GRAPHIC  OUTPUT  PILES.. WILL  ALTERNATELY 
MERGE  INTO  A  USER  CHOSEN  FILE  NAME. 

$  AMP  -  HISTOGRAMS  WITH  STATISTICS  FOR  SAMPLED  DATA  FILES 


$  CLR  PUTS  AN  ADM  TERMINAL  IN  ALPHA  MODE  THEN 
CLEARS  SCREEN  AFTER  GRAPHICS  SESSION 
$  CPF  -  REASSIGNS  GRAPHIC  OUTPUT  TO  THE  TERMINAL 
$  GRA  -  PUTS  AN  ADM  TERMINAL  IN  GRAPHICS  MODE 

THEN  CLEARS  SCREEN . 

$  CXP  -  MODIFICATION  OF  CEP  TO  WRITE  DATA  WHICH  WAS 
PREVIOUSLY  PLOTTED,  TO  A  FILE. 


$  BED  -  LIST  CONTENTS  OF  THE  SPS  PORTION  OP  ILS  HEADER 

$  HFL  -  APPEND  HEADER  TO  YOUR  DATA 

$  HPF  -  PROCESSES  GRAF.TMP  FILES  RECOGNIZING  FORM  FEED 
NOT  TO  BE  USED  WITH  OVERLAYS.  DEVELOPED  FOR 
USE  WITH  HSP. 

$  HSP  -  HISTOGRAM  OF  DATA  PEAKS  AND  ENERGY  CONTENT  OF 
DATA  FILES. 

$  ITN  -  INSERTS  TEST  NAME , COMMENTS  AND  ISPS  HEADER 
VALUES  INTO  ISPS  PORTION  OF  HEADER 
$  LAB  -  LABELS  FOR  PLOTS.  CREATES  A  GRAF.TMP  FILE  TO 
BE  MERGED  TO  PROVIDE  HORZ  _VERT  LABELING. 

$  MST  -  WRITES  STATISTICS  FROM  SAMPLED  DATA  FILES 
AS  FROM  ILS  COMMAND  "$  CST",  TO  THE  FILE 
"CSTAT.TMP"  WITH  LEGEND  AND  HEADINGS. 

SUBSEQUENT  USES  APPEND  DATA. 

$  OPP  -  ASSIGNS  GRAPHIC  OUTPUT  TO  A  FILE  NAMED  GRAF.TMP 
$  PPF  -  PROCESSES  AND  PRINTS  GRAF.TMP  FILES 

$  RR  -  LISTS  COMMANDS  WHICH  USE  RECORD  DATA 

$  RST  -  WRITES  STATISTICS  FROM  RECORD  DATA  PILES 

AS  FROM  ILS  COMMAND  "$  CST",  TO  THE  FILE 
"CSTAT.TMP"  WITH  LEGEND  AND  HEADINGS. 

SUBSEQUENT  USES  APPEND  DATA. 

$  RTS  -  RECORD  DATA  TO  SAMPLED  DATA  FOR  DISPLAY  PURPOSES 

$  STN  -  SEARCHES  A  DIRECTORY  FOR  WD  FILES  PROM  A 

GIVEN  TEST,  ALTERNATELY  USES  WD  NUMBERS  FROM 
A  FILE. 

$  TRM  -  TRM  M  WILL  PRINT  VERTICALLY  IF  Nl  IS  NEGATIVE 
$  WIN  -  APPLIES  FULL  OR  PARTIAL  HANNING  WINDOW 
TO  SAMPLED  DATA 

$  WTN  -  WRITES  ISPS  PORTION  OP  FILE  HEADER  TO 
SECONDARY  FILE 


$  XCP  -  COMPUTES  COMPLEX  CEPSTRUM  FROM  INPUT  WAVEFORM, 
(INCLUDES  EXPONENTIAL  WEIGHTING,  AND  BANDPASS 
MAPPING) 


Specialized  Interactive  "Recipes"  using  ILS  commands 

***************************************************************** 


***************************************************************** 


**  4t  **  4r  ****  4t  *********  4r  ***  4t  *************  4r  ****  4r  ****  4t  4t  **************  * 

XCPILOW.COH  -  an  ILS  "recipe”  that  implements  a  forward  complex 
cepstrumr  low-pass  filters  the  cepstrum/  and  finally  does  an 
inverse  complex  cepstrum. 

***************************************************************** 


$1  LOW-PASS  CEPSTRUM  AND  INVERSE  TRANSFORM 
$  ON  CONTROLY  THEN  GOTO  CLEANUP 


$  ON  ERROR  THEN  GOTO  CLEANUP 


$ 

IF 

PI  .EQS. 

m  a 

THEN 

INQUIRE 

$ 

IF 

P2  .EQS. 

m  m 

THEN 

INQUIRE 

$ 

IF 

P3  .EQS. 

m  m 

THEN 

INQUIRE 

$ 

IF 

P4  .EQS. 

m  m 

THEN 

INQUIRE 

$ 

IF 

P5  .EQS. 

m  m 

THEN 

INQUIRE 

$  CTX  'P5' 

$  FIL  SO 

CZ9996,,1 

CZ9997,,1 

$  FIL  'PI' 

$  FIL  S9996 
$  XCP  Z'P3' , •P4' 

$  FIL  9996 
$  FIL  S9997 
$  XCP  IS 

$  RENAME  WD9997.  WD'P2'. 
$CLEANUP: 

$  dele/NOCONFIRM  WD9996.;1 


PI  "INPUT  FILE  NUMBER?" 

P2  "OUTPUT  FILE  NUMBER?" 

P3  "STARTING  CHANNEL  NO.?" 
P4  "ENDING  CHANNEL  NO.?" 

P5  "CONTEXT?" 


***************************************************************** 
XCPIHIGH.COM  -  an  ILS  "recipe"  that  implements  a  forward  complex 
cepstrum,  high-pass  filters  the  cepstrum,  and  finally  does  an 
inverse  complex  cepstrum. 

***************************************************************** 

$1  HIGH-PASS  CEPSTRUM  AND  INVERSE  TRANSFORM 
$  ON  CONTROLY  THEN  GOTO  CLEANUP 
$  ON  ERROR  THEN  GOTO  CLEANUP 


$ 

IP 

PI  .EQS. 

m  n 

THEN 

INQUIRE 

PI 

"INPUT  FILE  NUMBER?" 

§ 

IP 

P2  .EQS. 

m  m 

THEN 

INQUIRE 

P2 

"OUTPUT  PILE 

NUMBER?" 

$ 

IF 

P3  .EQS. 

m  a 

THEN 

INQUIRE 

P3 

"U.P.  OP  1st 

BAND?" 

$ 

IP 

P4  .EQS. 

a  a 

THEN 

INQUIRE 

P4 

"L.B.  OF  2nd 

BAND?" 

$ 

IF 

P5  .EQS. 

a  a 

THEN 

INQUIRE 

P5 

"CONTEXT?" 

$  CTX  'PS' 
$  FIL  SO 
CZ9996,,1 
CZ9997,,1 
CZ9998,,1 


$  FIL  'PI' 

$  FIL  S9996 


$  FIL  9996 
$  FIL  S9997 
$  XCP  Z'P4' ,'P5' 

$  FIL  9997 
$  FIL  S9998 
$  XCP  IS 

$  RENAME  WD9998.  WD'P2'. 
$CLEANUP: 

$  dele/NOCONFIRM  WD9996.;1 
$  dele/NOCONFIRM  WD9997.;1 


XCPIPLOW.COM  -  an  ILS  "recipe”  that  implements  bandpass  mapping 
for  bandpass  limited  signals,  followed  by  a  forward  complex 
cepstrum,  a  low-pass  filtering  operation  on  the  cepstrum,  and 
finally  an  inverse  complex  cepstrum. 

It**************************************************************** 


$1  LOW-PASS  CEPSTRUM  AMD  INVERSE  TRANSFORM 
$1  MODIFIED  FOR  BAND-PASS  SYSTEM 
$  ON  CONTROLY  THEN  GOTO  CLEANUP 


$  ON  ERROR  THEN  GOTO  CLEANUP 


$ 

IF 

PI  .BQS. 

II  • 

THEN 

INQUIRE 

$ 

IF 

P2  .EQS. 

m  m 

THEN 

INQUIRE 

$ 

IF 

P3  .EQS. 

m  m 

THEM 

INQUIRE 

$ 

IP 

P4  .EQS. 

m  M 

THEN 

INQUIRE 

$ 

IF 

P5  .EQS. 

«  n 

THEN 

INQUIRE 

$  CTX  'P5' 

$  FIL  SO 

CZ9996,,1 

CZ9997,,1 

$  FIL  'PI' 

$  FIL  S9996 
$  XCP  PZ'P3' ,'P4' 

$  FIL  9996 
S  FIL  S9997 
$  XCP  ISP 

$  RENAME  WD9997.  WD'P2'. 
$CLEANUP; 

$  dele/NOCONFIRM  WD9996.;1 


PI  "INPUT  FILE  NUMBER?" 

P2  "OUTPUT  PILE  NUMBER?" 

P3  "STARTING  CHANNEL  NO.?" 
P4  "ENDING  CHANNEL  NO.?" 

P5  "CONTEXT?" 


***************************************************************** 
XCPIPHIGH.COM  -  an  ILS  "recipe"  that  implements  bandpass  mapping 
for  bandpass  limited  signals,  followed  by  a  forward  complex 
cepstrum,  a  high-pass  filtering  operation  on  the  cepstrum,  and 
finally  an  inverse  complex  cepstrum. 
***************************************************************** 


$1  HIGH-PASS  CEPSTRUM  AND  INVERSE  TRANSFORM 
$1  MODIFIED  FOR  BAND-PASS  SYSTEM 
$  ON  CONTROLY  THEN  GOTO  CLEANUP 


$  ON  ERROR  THEN  GOTO  CLEANUP 

$  IF  PI  .EQS.  ""  THEN  INQUIRE  Pi  "INPUT  FILE  NUMBER?" 

$  IF  P2  .EQS.  ""  THEN  INQUIRE  P2  "OUTPUT  FILE  NUMBER?" 

$  IF  P3  .EQS.  ""  THEN  INQUIRE  P3  "U.P.  OF  1st  BAND?" 

$  IF  P4  .EQS.  ""  THEN  INQUIRE  P4  "L.B.  OF  2nd  BAND?" 

$  IF  P5  .EQS.  ""  THEN  INQUIRE  P5  "CONTEXT?" 

$  CTX  'P5' 

$  FIL  SO 
CZ9996,,1 
CZ9997,,1 
CZ9998,,1 

$  FIL  'PI' 

$  FIL  S9996 
$  XCP  PZl, 'P3' 

$  FIL  9996 
$  FIL  S9997 
$  XCP  PZ'P4' , 'P5' 

$  FIL  9997 
$  FIL  S9998 
$  XCP  ISP 

$  RENAME  WD9998.  WD'P2'. 

$CLEANUP : 

$  dele/NOCONFIRM  WD9996.;1 
$  dele/NOCONFIRM  WD9997.;1 


************  It**************  *ii********it****iiit****it***ii**iiti*ii****** 

FDECON.COM  -an  ILS  "recipe"  that  implements  Fourier  Deconvolution 

***************************************************************** 


$1  FOURIER  DECONVOLUTION 
$  ON  CONTROLY  THEN  GOTO  CLEANUP 
$  ON  ERROR  THEN  GOTO  CLEANUP 

$  IF  PI  .EQS.  ""  THEN  INQUIRE  Pi  "INPUT  A-FILE  NUMBER?" 

$  IF  P2  .EQS.  ""  THEN  INQUIRE  P2  "INPUT  B-FILE  NUMBER?" 

$  IF  P3  .EQS.  ""  THEN  INQUIRE  P3  "OUTPUT  FILE  NUMBER?" 

$  FIL  S9993 
$  OPN  S5 
$  FIL  'PI' 

$  SRE  1,1 
$  FIL  'P2' 

$  FIL  S9994 
$  SRE  1,1 
$  FIL  9993 
$  FIL  S9995 
$  FFT 

$  FIL  9994 
$  FIL  S9996 
$  FFT 

$  FIL  9995 
$  FIL  B9996 
$  FIL  S9997 
$  BOP  D 


$  FIL  S'P3' 

$  OFN  SI 
$  PFT  I 
$CLEANUP: 

$  dele/NOCONFIRH  WD9993.;1 
$  dele/NOCONFIRM  WD9994.;1 
$  dele/NOCONFIRM  WD9995.;1 
$  dele/NOCONFIRM  WD9996.;1 
$  dele/NOCONFIRM  WD9997.;1 

*ik*****4t4t*****4t*********4t4t*44t*********4t****4t****************4!*itilr* 

FDECONP.COM  -  sn  ILS  "recipe"  that  implements  Fourier 

Deconvolution  on  bandpass  mapped  waveforms. 
***************************************************************** 

$1  FOURIER  DECONVOLUTION  WITHIN  A  FREQUENCY  BAND 
$  IF  PI  .EQS.  ""  THEN  INQUIRE  Pi  "INPUT  A-FILE  NUMBER?" 

$  IF  P2  .EQS.  ""  THEN  INQUIRE  P2  "INPUT  B-FILE  NUMBER?" 

$  IF  P3  .EQS.  ""  THEN  INQUIRE  P3  "OUTPUT  FILE  NUMBER?" 

$  IF  P4  .EQS.  ""  THEN  INQUIRE  P4  "FREQUENCY  LOWER  LIMIT?" 

$  IF  P5  .EQS.  ""  THEN  INQUIRE  P5  "FREQUENCY  UPPER  LIMIT?" 

$  ON  CONTROLY  THEN  GOTO  DELETECOMMAND 
$  OPEN/WRITE  OUT  FDE.COM 

$  WRITE  OUT  "$  ON  CONTROLY  THEN  GOTO  CLEANUP" 

$  WRITE  OUT  "$  ON  ERROR  THEN  GOTO  CLEANUP" 

$  WRITE  OUT  "$  FIL  S9993" 

$  WRITE  OUT  "$  OPN  S5" 

$  WRITE  OUT  "$  FIL  "PI'" 

$  WRITE  OUT  "$  SRE  1,1" 

$  WRITE  OUT  "$  FIL  "P2'" 

$  WRITE  OUT  "$  FIL  S9994" 

$  WRITE  OUT  "$  SRE  1,1" 

$  WRITE  OUT  "$  FIL  9993" 

$  WRITE  OUT  "$  FIL  S9995" 

5  WRITE  OUT  "$  FFT" 

$  WRITE  OUT  "$  FIL  9994" 

$  WRITE  OUT  "$  FIL  S9996" 

$  WRITE  OUT  "S  FFT" 

$  WRITE  OUT  "$  FIL  9995" 

$  WRITE  OUT  •$  FIL  B9996" 

$  WRITE  OUT  "$  FIL  S9997" 

$  WRITE  OUT  "$  BOP  DP" 

$  WRITE  OUT  "  "P4'.,"P5'." 

$  WRITE  OUT  "$  FIL  9997" 

$  WRITE  OUT  "$  FIL  S"P3'" 

$  WRITE  OUT  "$  OPN  Si" 

$  WRITE  OUT  "$  FFT  I" 

$  WRITE  OUT  "$CLEANUP:" 

$  WRITE  OUT  "$  dele/NOCONFIRM  WD9993.;1" 

$  WRITE  OUT  "$  dele/NOCONFIRM  WD9994.;1" 

$  WRITE  OUT  "$  dele/NOCONFIRM  WD9995.;1" 

$  WRITE  OUT  "$  dele/NOCONFIRM  WD9996.,*!" 

$  WRITE  OUT  "$  dele/NOCONFIRM  WD9997.}!" 

$  CLOSE  OUT 
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$  @FDE.COM 
$DELETECOHHAND : 

$  DELETE/NOCONFIRM  PDE.COM; 


***************************************************************** 
AVERG.COM  -  an  ILS  "recipe*  that  implements  an  averaging 
operation  ona  series  of  candidate  waveforms. 
***************************************************************** 


$  ON  CONTROLY  THEN  GOTO  CLEANUP 
$  ON  ERROR  THEN  GOTO  CLEANUP 

$  IF  PI  .EQS.  ""  THEN  INQUIRE  Pi  "INPUT  FILE?" 

$!  INPUT  FILE  CONTAINS  WD  NUMBERS,  ONE  PER  LINE,  OF  FILES  TO  BE 
AVERAGED 

$  IP  P2  .EQS.  ""  THEN  INQUIRE  P2  "OUTPUT  FILE  NUMBER?" 

$  OPEN/READ  IN  'Pi' 

$  COUNT* 0 
$  FIL  S9995 
$  OPN  S3 
SLOOP: 

$  READ/ENDOFFILE=THATSALL  IN  NUM 
$  C0UNT=C0UNT+1 
$  FIL  'NUM' 

$  SRE  1,1 
$  GOTO  LOOP 
STHATSALL: 

$  FIL  9995 
$  FIL  S9996 
$  AVG  1, 'COUNT' 

$  FIL  9996 
$  FIL  S9997 
$  TRE  ' COUNT ',1 
$  FIL  9997 
$  FIL  S’P2' 

$  RTS  1,1 
SCLEANUP: 

$  dele/NOCONFIRM  WD9995.;1 
$  dele/NOCONFIRM  WD9996.;1 
$  dele/NOCONFIRM  WD9997.;1 
$  CLOSE  IN 


***************************************************************** 

Special  Purpose  ILS  Software  Modules  created  for  this  Study 

***************************************************************** 


***************************************************************** 

XCP.POR 

***************************************************************** 
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n n o o n  n o n n  o o n n  o o r> n o  n n o  n o o o n o n o o n o o  o  oooo 


INTERACTIVE  LABORATORY  SYSTEM 
ILS  COMMAND  PROGRAM  **  XCP  ** 
PROGRAM  XCP 
IMPLICIT  INTEGER  (I-N) 

START  OF  DOCUMENTATION 

SPECIAL  VERSION  TO  USE  ENTIRE  FILE 


COMMAND  FORMAT; 

XCP  [Z,I,P] [S]N1,N2 
ALPHABETIC  ARGUMENTS: 

Z  -  ZERO  CEPSTRUM  RECORD  FROM  Nl  TO  N2 

I  -  INVERSE  CEPSTRUM 

F  -  CREATE  TEST  FILE 

S  -  USED  WITH  [I],  OPTIONAL  PHASE  SHIFT 

P  -  PREPROCESS  WITH  FREQUENCY  LIMIT 

NUMERIC  ARGUMENTS; 

Nl  -  STARTING  FRAME 
N2  -  NUMBER  OF  FRAMES 
WITH  OPTION  [ZJ 
Nl  -  FIRST  POINT  TO  ZERO 
N2  -  LAST  POINT  TO  ZERO 

N3  -  ASK  FOR  ALPHA 
WITH  OPTION  [P] 

N4  -  LOWER  FREQUENCY 
N5  -  UPPER  FREQUENCY 

END  OF  DOCUMENTATION 

COMMON 

PI , TWOPI , THLINC , THLCON , NFFT , NDUM , NN , L , H , HI , DVTMN2 , AMULT 
COMMON  /CLBF/  INSFLG ,LCLBF , ICLBF ( 40) 

COMMON  /ILSA/  NBCW,NCWBK,NDPBK,NDPF,NBDP,NCWFH, 

1  KBU , KBUIN , LPU , LUGI , LUGO , NSC , CWSC , 

2  NBA2D,MIDA2D,ICTIM{4) ,ICDAT(6) 

COMMON  /ILSB/  IA( 4) ,Nl ,N2 ,N3 ,N4 ,N5 ,N6 ,N7 ,N8 ,N9 ,N10 ,Nll ,N12 
COMMON  /ILSC/  IASAV(4) ,NlSAV,N2SAV,N3SAVrN4SAV,N5SAV,N6SAV, 
1  N7  SAV , N8SAV , N9SAV , Nl 0SAV , Nl ISAV , Nl 2SAV 

COMMON  /ILSE/ 

I F  LPA (16), LENPA , LF I LPA , NFLP A , I DKPA , I DKDP A , I DP A 
COMMON  /ILSF/ 

IFLSA ( 16 )  r  LENSA , LFILSA , NFLSA , IDKSA , IDKDSA , IDSA 

COMMON  /ILSH/  FS ,M ,MPl ,M02 ,N,NSPBK ,NSHFT , ICON, ISP , IHAM , LAN 
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,  ■ : 


oo  no  on  on  non  non  n  non 


COMMON  /ILSI/  LRH,IEX,ISTAN,NAN,IDF(5) ,NP,IVL,IC 
COMMON  /ILSJ/  RC(30) ,A(30) ,R(30) ,F(10) ,B(10) 

C. .  • 

COMMON  /FLCl/  IFLl ( 16) ,LENl ,LFILl ,NFLl , IDKl , IDKDl , IDl 
COMMON  /FLCM/ 

IFLC(16) ,LENC,IFLALF(4) ,LENALF,IFLX,IFLV,LENCOM 

C . . . 

DIMENSION  IY(16384) ,IR(16896) ,IS(16896) ,IX(16384) 
DIMENSION  Y(16384) ,X(16384) ,DUM(16384) 

DIMENSION  IHR(128) ,IHW(128) 

DIMENSION  I0PT(5) ,ICHSTR(5) 

EQUIVALENCE  (IHR(70) ,SCALR) ,(IHW(70) ,SCALW) 
EQUIVALENCE  ( IHR(74) ,ACEP)  ,( IHW(74) , ALPHA) 


,2HF  / 


C  •  •  • 


DATA  SCLl/1.0/ 

DATA  11,12,13,14,15,111,119,130/1,2,3,4,5,11,19,30/ 

DATA  ICHSTR(l) ,ICHSTR(2) ,ICHSTR(3) ,ICHSTR(4)/2HZ  ,2HI  ,2HS 

DATA  ICHSTR(5)/2HP  / 

DATA  IHW(75) ,IHW(76)/0,0/ 

DATA  IFLG , I IBW , I IBR , IOVERFLO/4 *  0/ 

DATA  LIX,LIY,LIR,LIS/2*16384, 2*16896/ 

DATA  ICHA,ICHNP,ICHS,ICHR/2HA  ,2HNP,2HS  ,2HR  / 

DATA  SCL0/0.0/ 

CALL  RCOMM 


INITIALIZATION  DATA 

CALL  AFARG ( lA , I 1 , I 4 , ICHSTR , lOPT ,15) 

CHECK  OPTIONS 

IF(IOPT(l)  .EQ.DTHEN 

ZERO  PORTION  OF  CEPSTRUM 
IMODE=l 

ELSE  IF  (IOPT(2)  .EQ.DTHEN 

INVERSE  CEPSTRUM 
IMODE=2 

ELSE  IF  (IOPT(4)  .EQ.DTHEN 

CREATE  TEST  FILE 
IMODE=4 

ELSE 

FORWARD  CEPSTRUM  OR  PREPROCESS 
IMODE=0 

END  IF 


/•  ••  .*►  '•* 


INVERSE  CEPSTRUM  WITH  PHASE  SHIFT 
IF(IMODE.EQ.2.AND.IOPT(3) .EQ.l) IM0DE=3 

SET  FRAMES  TO  DEFAULT  TO  WHOLE  FILE 

DO  10  1=1,4 

IF(IOPT(I) .NE.0)GOTO  11 
CONTINUE 

FIND  CEPSTRUM . SET  UP  NUMBER  OF  POINTS 

IP(N1.EQ.0.AND.N2.EQ.0)THEN 

N1A=0 

N2A=1 

ELSE 

N1A=N1 

N2A=N2 

END  IF 
GOTO  12 

OPTIONS  [I]  AND  [Z]  USE  THE  ENTIRE  FILE 
IF(IOPT(4)  .EQ.DTHEN 
NPTS=2048 
LFILSA=2048/NDPBK 
LFRAM=NDPF 
NFRAM=2048/NDPF 
GOTO  260 

END  IF 

N1A=0 

N2A=1 

CALL  CHKFL ( NFLPA , IDKPA , IFLPA , LENPA , LFILPA , IHR , ICHS , lERR) 
IF(IERR.NE.0)  GO  TO  310 

GET  SAMPLING  FREQUENCY 

ISF=IHR(62) 

IPWR=IHR(61) 

CALL  SFCNV(FS,ISF,IPWR,I1) 

CALL  ANCHK ( NSC , N1 A , N2 A , N1 SAV , N2 S AV , NSCA , I STAN , 

1  NAN, NDPF,NSDBK, LFILPA) 

ISTFR=N1SAV 

NPR=N2SAV 

GET  NUMBER  OF  POINTS . TRUNCATE  TO  A  POWER  OF  2 

CALL  FTOP ( N1 SAV , N2 SAV , NSC , NDPF , STPT , PTS , ENDPT ) 
NPTS=IFIX(PTS+0.5) 

INDX=16384 

IF(NPTS.GE.INDX)THEN 
NPTS=INDX 
GOTO  30 


ELSE 


INDX=INDX/2 


END  IF 
GOTO  20 


CONVERT  BACK  TO  FRAMES 
PTS=FLOAT{NPTS) 

CALL  PTOF ( STPT , PTS , NDPF , NSC , ISTFR , NFR) 

NPTX=NPTS 
NPTY»NPTS 
IST»ISTFR 
LST»ISTFR+NFR-1 
SET  SCALE 

IF ( INODE. EQ.l) THEN 
SCALE=1 . 0 

ELSE  IF(IHODE.EQ.0)THEN 
SCALE^l . 0 

ELSE  IF ( INODE . EQ . 2 .OR . INODE . EQ . 3) THEN 
SCALEs^SCALR 

END  IF 

CALL  GETD ( NSC , I STFR , NPTS , NPTY , IR , LIR , I Y , I IBR , 

1  IFLPA,LENPA,LFILPA) 

IF(IIBR.LT.0)GOTO  260 

CALL  MI2R( I Y,Y, NPTS, SCALE) 

CONTINUE 

ALL  POINTS  READ 

LFILSA«LFILPA 

IPLGO=2 

ITYPE=-1 

CALL  CHKFL ( NFLSA , IDKSA , IFLSA , LENSA , LFILSA , IHW , ITYPE , IFLGO ) 
ITYPE=ICHNP 

CALL  SETUP ( IFLSA , LENSA , LFILSA , NFRAM , LFRAM , NCWFH , I HW , ITYPE ) 
IP(NFRAM.EQ.-l)GOTO  310 


IHW(62)  *  IHR(62) 

IHW(61)  »  IHR(61) 

IHW(63)  «  -32000 
IP(IOPT(5)  .EQ.DTHEN 

IF (INODE. EQ.0) THEN 
Fl*FLOAT(N4) 
F2=FLOAT(N5) 
IHW(75)=N4 
IHW(76)=N5 

ELSE 

Fl=FLOAT(IHR(75) ) 
F2=FLOAT(IHR(76) ) 
IHW(75)»IHR(75) 
IHW(76)=IHR(76) 

END  IF 

ELSE 

IHW(75)=0 

IHW(76)»0 


P 


r 

r 
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END  IF 

IF ( INODE . EQ . 1 ) SCALWaSCALR 
C 

NCEP=NPTS 

AMULT=ACEP 

0  •  •  • 

C...  HAP  FORWARD  TRANSFORM 

C . . . 

IF(IOPT(5) . EQ.l. AND. INODE. EQ.0) THEN 

CALL  CMAP(NCEP,Y,X,FS, FI, F2, INODE) 

C  STOP 

CALL  MR2R(X,Y,NCEP,I1) 

END  IF 

CALL  CCPC(NCEP,Y,X, INODE) 

C... 

C...  MAP  INVERSE??? 

C . . . 

IF(IOPT(5) .EQ.l.AND. (IM0DE.EQ.2.0R.IM0DE.EQ.3) )THEN 
CALL  CMAP(NCEP,X,Y,PS,Pl,F2,IMODE) 

CALL  MR2R(Y,X,NCEP,I1) 

END  IF 

C . . . 

C...  X  IS  THE  OUTPUT  VARIABLE 

C . . . 

C...  SCALE  OUTPUT  VARIABLE 

C  . . 

IF (INODE. EQ.l) GOTO  275 
CALL  FPPIC(X,Il,NPTS,JLOC) 

SCALE=ABS( 32767 ./X(JLOC) ) 

SCALW=1 ./SCALE 

C...  NOW  THAT  OUTPUT  SCALE  IS  KNOWN.. WRITE  HEADER 

C...  SCALE  IS  EQUIVALENCED  TO  HEADER  LOCATION  70 

C...  WHEN  THE  ZERO  OPTION  IS  USED,  NO  SCALING  IS  DONE 

C...  ALSO  SET  ALPHA=AMULT  TO  STORE  IN  HEADER 

C .  . . 

275  ALPHA*AMULT 

CALL  WHEAD(IHW,IFLSA,LENSA) 

C . . . 

C...  CONVERT  TO  INTEGER. .SCALE  FOR  MAXIMUM  ACCURACY 

C . . . 

CALL  MR2I(X, IX, NPTS, SCALE) 

IFRAM=1 

CALL  WRITD ( NSC , IFRAM , I SDB , NPTS ,IS,LIS,IX,IIBW, 

1  IFLSA,LENSA,LFILSA) 

280  CONTINUE 

C...  DUMP  BUFFER 

IIBW»-IIBW 

CALL  WRITD (NSC, IFRAM, ISDB, NPTS, IS, LIS, IX, IIBW, 

1  IFLSA,LENSA,LFILSA) 

100  CONTINUE 

999  CONTINUE 

310  CALL  WCOMH 

CALL  EXILS 
END 


L 
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List  of  Fortran  Modules  needed  to  compile  XCP.FOR 

********4t4t*4r*******4t******4t***4t********  **********  ******4^********* 


AMODSQ 

FFA 

FFS 

FFT 

FFT842 

MRIDF 

ORDl 

ORD2 

PHADVT 

PHAUNW 

PHCHCK 

PPVPHA 

R2TR 

R2TX 

R4SYN 

R4TR 

R4TX 

R8SYN 

R8TR 

R8TX 

RP 

SPCVAL 


contained  in  CCEPS.FOR  and  CCXTRA.FOR 


SHIFTF 


***************************************************************** 


CMAP.FOR 

***************************************************************** 


SUBROUTINE  CMAP ( NSAMP ,Y,B,FS,Fl,P2, MODE ) 

C.  .  . 

C...  MAPS  FILTERED  DATA  WITH  A  SELECTED  FREQUENC'”  RANGE 

C...  INTO  ENTIRE  SPECTRUM  FOR  FORWARD  CEPSTRUM  OR  ENTIRE 

C...  FREQUENCY  RANGE  INTO  SELECTED  BAND  FOR  INVERSE  CEPSTRUM 

C. .  . 

C...  ARGUMENTS: 

C . . .  NSAMP 

C.  .  .  Y 

C  •  •  •  B 

C...  PS 

C...  PI 

C...  F2 

C . . .  MODE 

C... 

^* *************************************************************** 


-  NUMBER  OP  SAMPLES 

-  INPUT  VARIABLE  TIME  SERIES 

-  OUTPUT  VARIABLE  TIME  SERIES 

-  SAMPLING  FREQUENCY 

-  LOWER  FREQUENCY  BOUND 

-  UPPER  FREQUENCY  BOUND 

-  I MODE  PROM  CEPSTRUM 


C... 
C... 
C... 
C... 
C..  . 


VARIABLE  LIST 

DELF  FREQUENCY  INCREMENT 

DELT  TINE  INCREMENT 

DINT  INTERPOLATING  FRACTION 

PINT  REAL  LOCATION  OF  LOWER  FREQUENCY  BOUND 


C...  FMAX  MAXIMUM  FREQUENCY 

C...  ILOC  INTEGER  LOCATION  OF  INTERPOLATED  FREQUENCY 

C...  INTFl  INTEGER  LOCATION  OF  LOWER  FREQUENCY  BOUND 

C...  INTF2  INTEGER  LOCATION  OF  UPPER  FREQUENCY  BOUND 

C...  lOFSET  ADDER  TO  LOCATE  FREQUENCY  RANGE 

C...  IPOS  LOCATION  OF  REAL  PART  OF  LOWER  INTERPOLATION  PAIR 

C...  NINT  NUMBER  OF  FREQUENCIES  TO  MAP{F2-Fl),  DOES  NOT 

C  INCLUDE  FI 

C...  R1  REAL  PART  OP  LOWER  INTERPOLATION  PAIR 

C..  R2  REAL  PART  OF  UPPER  INTERPOLATION  PAIR 

C...  RATIO  RATIO  OF  MAPPED  FREQUENCIES  TO  TOTAL  FREQUENCIES 

C...  VLOC  LOCATION  OF  INTERPOLATED  FREQUENCY 

C...  XI  IMAGINARY  PART  OP  LOWER  INTERPOLATION  PAIR 

C...  X2  IMAGINARY  PART  OF  UPPER  INTERPOLATION  PAIR 

(2* ******************************* ******************************** 
DIMENSION  Y(l) ,B(NSAMP+2) 

delt=i./fs 

DELF=FS/NSAMP 

NFREQ=NSAMP/2. 

FMAX=FS/2. 

C.  .  . 

C...  WHERE  IN  THE  FREQUENCY  ARRAY  IS  THE  LOWER  FREQUENCY  BOUND 

C  .  • 

FINT=F1/DELF 

INTF1=IFIX(FINT) 

IF (FLOAT (INTFl) .LT.FINT) INTFl*INTFl+l 
INTF2=IFIX ( F2/DELF) 

NINT=INTF2-INTF1 

IOFSET=2*INTFl+l 

C  • 

C*. .  CONVERT  TO  FREQUENCY  DOMAIN 

C  •  •  • 

c 

C...  SUBROUTINE  FFA  REPLACES  THE  REAL  VECTOR  B(K) , 

( K  =  1  f  2  f  ,  m  •  fN)  f 

C...  WITH  ITS  FINITE  DISCRETE  FOURIER  TRANSFORM.  THE  DC  TERM  IS 
C...  RETURNED  IN  LOCATION  B(l)  WITH  B(2)  SET  TO  0.  THEREAFTER, 
THE 

C...  JTH  HARMONIC  IS  RETURNED  AS  A  COMPLEX  NUMBER  STORED  AS 

C...  B(2*J+1)  +  I  B(2*J+2).  NOTE  THAT  THE  N/2  HARMONIC  IS 

RETURNED 

C...  IN  B(N+1)  WITH  B(N+2)  SET  TO  0.  HENCE,  B  MUST  BE 
DIMENSIONED 
C. . .  TO  SIZE  N+2. 

C...  SUBROUTINE  IS  CALLED  AS  FFA  (B,N)  WHERE  N=2**M  AND  B  IS  AN 
C...  N  TERM  REAL  ARRAY. 

C  .  . 

CALL  PFA(Y,NSAMP) 

C 

C...  Y  NOW  CONTAINS  FREQUENCY  VALUES 

C . . .  MAP  SELECTED  RANGE  INTO  ARRAY  B 

C .  .  . 

IF ( MODE. EQ.0) THEN 

B(l) =Y(IOFSET) 

B(2) =Y(IOPSET+l) 


RATIO=FLOAT ( MINT) /FLOAT ( NFREQ) 
TYPE  *,DELT,DELF, NFREQ, FMAX, RATIO 
TYPE  *,INTF1,INTF2,NINT,I0FSET 
DO  100  KJ=1, NFREQ 

VLOC=FLOAT{KJ) *RATIO 

ILOC»IFIX(VLOC) 

DINT*VLOC-ILOC 

IPOS=2*II.OC+IOFSET 

Rl=Y(IPOS) 

R2=Y(IPOS+2) 

Xl=Y(IPOS+l) 

X2»Y(IPOS+3) 

B(2*KJ+1) -Rl+DINT* (R2-R1) 
B( 2*KJ+2) =X1+DINT* (X2-X1) 


CONTINUE 

FORMAT( (I8,4(4F7.0,2X) )/) 
SHIFT  PHASE 


CALL  SHIFTF(B,NSAMP) 

WRITE ( 3 , 9997 ) 

FORMAT (8X, 4 (5X, 'INPUT •9X, 'SHIFTED' ,4X) ) 

WRITE(3,9998) 

FORMATC  FREQ  ',4(2('  REAL  IMAG  '))) 
WRITE(3,9999) ( J, (Y(2*J+1) ,Y(2*J+2) ,B(2*J+1) ,B(2*J+2) 

1  Y(2*J+3) ,Y(2*J+4) ,B(2*J+3) ,B(2*J+4) , 

2  Y(2*J+5) ,Y(2*J+6) ,B(2*J+5) ,B(2*J+6) , 

3  Y(2*J+7) ,Y(2*J+8) ,B(2*J+7) ,B(2*J+8) ) 

4  J=0,NFREQ-3,4) 

ELSE 

RATIO=FLOAT ( NFREQ) /FLOAT ( NINT) 

B(1)=0. 

B(2)=0. 

DO  200  KJ=1, NFREQ 

IF(KJ.LT.INTF1.0R.KJ.GT.INTF2)THEN 

B(2*KJ+1)=0. 

B(2*KJ+2) =0. 

PT  QP 

INDX=KJ-INTF1 

VLOC=RATIO*FLOAT( INDX) 

ILOC=IFIX(VLOC) 

DINT=VLOC-ILOC 

IPOS=2*ILOC+l 

Rl*Y(IPOS) 

R2=Y(IP0S+2) 

Xl*Y(IPOS+l) 

X2*Y(IPOS+3) 

B(2*KJ+1) =R1+DINT*(R2-R1) 
B(2*KJ+2) =X1+DINT*(X2-X1) 

END  IF 

CONTINUE 


END  IF 


oooonnnnn 


w  ««***  4r  ********  «r  *******  4r  *********  A  ********  4r  *****************  Hr  **  * 

...  SUBROUTINE  FFS  SYNTHESIZES  THE  REAL  VECTOR  B(K) ,  WHERE 
...  K»1,2,...,N.  THE  INITIAL  FOURIER  COEFFICIENTS  ARE  PLACED  IN 

...  THE  B  ARRAY  OF  SIZE  N'l-2.  THE  DC  TERN  IS  IN  B(l)  WITH 

...  B(2)  EQUAL  TO  0. 

...  THE  JTH  HARMONIC  IS  STORED  AS  B(2*J+1)  +  I  B(2*J+2). 

...  THE  N/2  HARMONIC  IS  IN  B(N+1)  WITH  B(N+2)  EQUAL  TO  0. 

...  THE  SUBROUTINE  IS  CALLED  AS  FFS{B,N)  WHERE  N»2**N  AND 
...  B  IS  THE  N  TERM  REAL  ARRAY  DISCUSSED  ABOVE. 

CALL  FFS(B,NSAMP) 

C.. 

RETURN 

END 

It**************************************************************** 

CCEPS.FOR 

***************************************************************** 


C  SUBROUTINE:  CCEPS 

C  SUBROUTINE  TO  COMPUTE  THE  COMPLEX  CEPSTRUM  OF  A  SEQUENCE  X(N) 


SUBROUTINE  CCEPS ( NX , X , I SNX , I SFX , I SSUC , CX , AUX ) 

DIMENSION  X(l) ,CX(1) ,AUX(1) 

COMMON  PI , TWOPI , THLINC , THLCON , NFFT , NPTS , N , L , H , Hi , DVTMN2 

LOGICAL  issue 

NPTS=NPFT/2 

N=12 

L=2**N 

H=FLOAT(L) *FLOAT(NFFT) 

H1=PI/H 

ISSUC=.TRUE. 

ISNX=1 


DO  10  I»1,NX 
CX(I)=X(I) 

AUX(I)*FLOAT(I“l)  *X(I) 
10  CONTINUE 
INITL=NX+1 
IEND*NFFT+2 
DO  20  I=INITL,IEND 
CX(I)*0.0 
AUX(I)»0.0 
20  CONTINUE 
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CALL  PPA(CX,NFFT) 

CALL  PFA(AUX,NPFT) 

IP(CX(1) .LT.0.0) ISNX=-1 

I0«-1 

DVTMN2=0.0 
IEND=NPTS+1 
DO  30  I=1,IEND 
10=10+2 
IE-IO+1 

AMAGSQ=AHODSQ(CX(IO) ,CX(IE)) 

PDVT=PHADVT(CX(IO) ,CX(IE) ,AUX(IO) ,AUX(IE) ,AMAGSQ) 

AUX ( 10) =AHAGSQ 
AUX(IE)»PDVT 
DVTMN2*DVTMN2+PDVT 
30  CONTINUE 

DVTMN2* ( 2 . *DVTMN2-AUX (2) -PDVT) /FLOAT (NPTS) 

PPDVT»A0X(2) 

PPHASE=0 . 0 

PPV=PPVPHA(CX(1) ,CX(2) ,ISNX) 

CX(1)=.5*AL0G(AUX(1) ) 

CX(2)=0.0 

10=1 

DO  50  I=2,IEND 
10=10+2 
IE=I0+1 
PDVT=AUX(IE) 

PPV=PPVPHA(CX(IO) ,CX(IE) rISNX) 

PHASE=PHAUNW ( X , NX , I SNX , I , PPHASE , PPDVT , PPV , PDVT , I SSUC ) 

IF (issue) GO  TO  40 
ISSUC=. FALSE. 

RETURN 

40  PPDVT=PDVT 

PPHASE=PHASE 

ex (10) = . 5*AL0G ( AUX ( 10) ) 

CX(IE)=PHASE 
50  CONTINUE 

I SPX=(ABS( PHASE/PI) +.1) 

IP (PHASE.lt. 0.0) ISFX=-ISFX 
H=PHASE/FLOAT ( NPTS ) 

IE=0 

DO  60  I=1,IEND 
IE-IE+2 

ex (IE) =CX(IE)-H*FL0AT(I-1) 

60  CONTINUE 


n  n 


CALL  FPS(CX,NPFT) 

RETURN 

END 


C 

C 


C  SUBROUTINE;  SPCVAL 

C  SUBROUTINE  TO  COMPUTRE  A  SPECTRAL  VALUE  AT  A  FREQUENCY 
C  FREQ(RAOIANS)  FOR  SEQUENCE  X(N)  AND  N*X(N) 

C - 


c 

SUBROUTINE  SPCVAL ( NX , X , FREQ , XR , XI , YR , YI ) 

DIMENSION  X(l) 

DOUBLE  PRECISION  U0 ,Ul ,U2 ,W0 ,W2,A,B,C,D,Al ,A2 ,SA0 ,CA0 


CA0-DBLE ( COS ( FREQ) ) 

SA0=DBLE( SIN (FREQ) ) 

A1»2.D+0*CA0 

U1S0.D+0 

U2=U1 

W1=U1 

W2=U1 

C 

DO  10  J=1,NX 
XJ=DBLE(X(J) ) 

U0*XJ+A1*U1-U2 

W0=DBLE ( FLOAT (J-l) ) *XJ+A1*W1-W2 

U2=U1 

U1=U0 

W2=W1 

W1=W0 

10  CONTINUE 
C 

A=U1-U2*CA0 

B=U2*SA0 

C=W1-W2*CA0 

D=W2*SA0 

A2=DBLE(FREQ*FLOAT(NX-l)  ) 
Ul=:DCOS(A2) 

U2*-DSIN(A2) 

XR*SNGL(U1*A-U2*B) 

XI=SNGL(U2*A+U1*B) 

YR=SNGL(U1*C-U2*D) 

YI»SNGL(U2*C+U1*D) 

RETURN 

END 

C 

C - 


o  o 


THE  UNWRAPPED  PHASE  ESTIMATE  IS  RETURNED  IN  PHAUNW. 


C 

FUNCTION  PHAUNW ( X , NX , I SNX , I , PPHASE , PPDVT , PPV , PDVT , I SCONS ) 
C 

DIMENSION  SDVT(17) ,SPPV{17) ,X(1) 

INTEGER  SINDEX(17) ,PINDEX,SP 
LOGICAL  ISCONS  FIRST 

COMMON  PI , TWOPI , THLINC , THLCON , NPPT , NPTS , N , L , H , Hi , DVTMN2 
C 

FIRST*. TRUE. 

PINDEX*! 

SP=1 

SPPV(SP)=PPV 

SDVT(SP)=PDVT 

SINDEX(SP)*L+1 

C 

GO  TO  40 
C 

10  PINDEX=SINDEX(SP) 

PPHASE=PHASE 
PPDVT*SDVT ( SP) 

SP*SP-1 
GO  TO  40 
C 

20  IF(  (SINDEX(SP)-PINDEX)  .GT.DGO  TO  30 
ISCONS*. FALSE. 

PHAUNW*0. 

RETURN 

C 

30  K*(SINDEX(SP)+PINDEX)/2 
C 

PREQ*TWOPI*( FLOAT (1-2) *FLOAT(L) +FLOAT ( K-1) ) /H 
CALL  SPCVAL ( NX , X , FREQ , XR , XI , YR , Y I ) 

C 

SP*SP+1 

SINDEX(SP)*K 

SPPV(SP) *PPVPHA(XR,XI,ISNX) 

XNAG*AMODSQ ( XR , XI ) 

SDVT ( SP ) * PH AD VT (XR,XI,YR,YI, XM AG ) 

C 

40  DELTA*Hl*FLOAT(SINDEX(SP) -PINDEX) 

PHAINC*DELTA* ( PPDVT+SDVT ( SP) ) 

C 

IP (ABS (PHAINC-DELTA*DVTMN2) .GT. THLINC) GO  TO  20 
C 

PHASE*PPHASE+FHAINC 

CALL  PHCHCK( PHASE, SPPV(SP) , I SCONS) 

IF ( .NOT. ISCONS) GO  TO  20 
C 

IP(ABS(PHASE-PPHASE)  .GT.PDGO  TO  20 


3 
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C  FUNCTION;  PPVPHA 

C  COMPUTE  THE  PRINCIPLE  VALUE  OF  THE  PHASE  OP  A  SPECTRAL  VALUE 
C - 


C 

FUNCTION  PPVPHA(XR,XI,ISNX) 

C 

IP(ISNX.EQ.1)PPVPHA=(ATAN2( (XI) ,(XR)  )  ) 
IPdSNX.EQ.  (-1)  )PPVPHA=(ATAN2(-(XI)  ,-{XR)  )  ) 
RETURN 
END 
C 

C - 


C  FUNCTION;  PHADVT 

C  COMPUTE  THE  PHASE  DERIVATIVE  OF  A  SPECTRAL  VALUE  OP  A  SEQUENCE 
X(N) 

C - 


c 

FUNCTION  PHADVT(XR,XI,YR,YI,XMAG) 
C 


PHADVT=-SNGL( (DBLE(XR) *DBLE (YR) +DBLE(XI) *DBLE ( YI) ) /DBLE (XMAG) ) 
RETURN 
END 
C 

C - 


C  FUNCTION;  AMODSQ 

C  COMPUTE  THE  SQUARE  OF  THE  MODULUS  OF  A  COMPLEX  NUMBER 

C - 


C 

FUNCTION  AMODSQ(ZR,ZI) 

C 

AMODSQ«SNGL(DBLE(ZR)  *DBLE  ( ZR) -t-DBLE  (  ZI )  *DBLE(ZI)  ) 
RETURN 
END 
C 

C - 


C  SUBROUTINE;  PHCHCK 

C  SUBROUTINE  TO  CHECK  CONSISTENCY  OF  A  PHASE  ESTIMATE 

C - 


Jl _ i 
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C 

SUBROUTINE  PHCHCK(PH,PV,ISCONS) 

C 

COMMON  PI , TWOPI , THLINC , THLCON , NFFT , NPTS , N , L , H , HI , DVTMN2 
LOGICAL  I SCONS 
C 

A0«{PH-PV) /TWOPI 
A1*PLOAT(IFIX(A0) ) *TWOPI+PV 
A2«A1+SIGN(TWOPI,A0) 

A3=ABS{A1-PH) 

A4sABS(A2-PH) 

C 

ISCONS». FALSE. 

I F ( A3 . GT . THLCON . AND , A4 . GT . THLCON ) RETURN 
ISCONS=.TRUE. 

C 

PH*A1 

IF(A3.GT.A4)PH=A2 

RETURN 

END 

C 


C  SUBROUTINE;  FFA 
C  FAST  FOURIER  ANALYSIS  SUBROUTINE 

C - 

c 

SUBROUTINE  FFA(B,  NFFT) 

C 

C  THIS  SUBROUTINE  REPLACES  THE  REAL  VECTOR  B(K) ,  (K=l ,2 , . . . ,N) , 

C  WITH  ITS  FINITE  DISCRETE  FOURIER  TRANSFORM,  THE  DC  TERM  IS 
C  RETURNED  IN  LOCATION  B(l)  WITH  B(2)  SET  TO  0.  THEREAFTER,  THE 
C  JTH  HARMONIC  IS  RETURNED  AS  A  COMPLEX  NUMBER  STORED  AS 
C  B(2*J+1)  +  I  B(2*J+2).  NOTE  THAT  THE  N/2  HARMONIC  IS  RETURNED 
C  IN  B(N+1)  WITH  B(N+2)  SET  TO  0.  HENCE,  B  MUST  BE  DIMENSIONED 
C  TO  SIZE  N+2. 

C  SUBROUTINE  IS  CALLED  AS  FFA  (B,N)  WHERE  N=2**M  AND  B  IS  AN 
C  N  TERM  REAL  ARRAY.  A  REAL-VALUED,  RADIX  8  ALGORITHM  IS  USED 
C  WITH  IN-PLACE  REORDERING  AND  THE  TRIG  FUNCTIONS  ARE  COMPUTED  AS 
C  NEEDED. 

C 

DIMENSION  B(2) 

COMMON  /CON/  PII,  P7,  P7TWO,  C22,  S22,  PI2 

IW  IS  A  MACHINE  DEPENDENT  WRITE  DEVICE  NUMBER 

IW  *  I1MACH(2) 

IW=6 
C 

PII  =  4.*ATAN(1.) 

PI8  »  PII/8. 

P7  »  l./SQRT(2.) 
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P7TW0  *  2.*P7 
C22  »  COS (PIS) 

S22  »  SIN(PI8) 

PI2  »  2.*PII 
N  =  1 

DO  10  1=1,15 
M  =  I 
N  =  N*2 

IF  (N.EQ.NPFT)  GO  TO  20 
10  CONTINUE 

WRITE  (IW,9999) 

9999  FORMAT  (30H  NFFT  NOT  A  POWER  OF  2  FOR  FFA) 
STOP 

20  CONTINUE 

N8POW  =  M/3 


DO  A  RADIX  2  OR  RADIX  4  ITERATION  FIRST  IF  ONE  IS  REQUIRED 

IF  (M--N8POW*3-l)  50,  40,  30 
30  NN  =  4 

INT  =  N/NN 

CALL  R4TR(INT,  B(l) ,  B(INT+1),  B(2*INT+1),  B(3*INT+1)) 

GO  TO  60 
40  NN  =  2 

INT  =  N/NN 

CALL  R2TR(INT,  B(l) ,  B(INT+1)) 

GO  TO  60 
50  NN  »  1 

PERFORM  RADIX  8  ITERATIONS 

60  IF  (N8POW)  90,  90,  70 
70  DO  80  IT=l,N8POW 
NN  =  NN*8 
INT  =  N/NN 

CALL  R8TR(INT,  NN,  B(l) ,  B(INT+1),  B(2*INT+1),  B(3*INT+1), 

*  B(4*INT+1) ,  B(5*INT+1),  B(6*INT+1),  B(7*INT+1),  B(l), 

*  B(IFT+1),  B(2*INT+1),  B(3*INT+1),  B(4*INT+1), 
B(5*INT+1) , 

*  B(6*INT+1),  B(7*INT+1)) 

80  CONTINUE 

PERFORM  IN-PLACE  REORDERING 

90  CALL  ORDKM,  B) 

CALL  ORD2(M,  B) 

T  =  B(2) 

B(2)  =  0. 

B(NPFT+1)  =  T 
B(NFFT+2)  =  0. 

DO  100  1=4, NFFT, 2 
B(I)  =  -B(I) 


o  o  o  n 


100  CONTINUE 
RETURN 
END 
C 

C - 

C  SUBROUTINE:  FFS 

C  FAST  FOURIER  SYNTHESIS  SUBROUTINE 
C  RADIX  8-4-2 

C - 

C 

SUBROUTINE  FFS(B,  NFFT) 

C 

C  THIS  SUBROUTINE  SYNTHESIZES  THE  REAL  VECTOR  B(K),  WHERE 
C  K=1,2,...,N.  THE  INITIAL  FOURIER  COEFFICIENTS  ARE  PLACED  IN 
C  THE  B  ARRAY  OF  SIZE  N+2.  THE  DC  TERM  IS  IN  B(l)  WITH 
C  B(2)  EQUAL  TO  0. 

C  THE  JTH  HARMONIC  IS  STORED  AS  B(2*J+1)  +  I  B(2*J+2) . 

C  THE  N/2  HARMONIC  IS  IN  B(N+1)  WITH  B(N+2)  EQUAL  TO  0. 

C  THE  SUBROUTINE  IS  CALLED  AS  FFS(B,N)  WHERE  N=2**M  AND 
C  B  IS  THE  N  TERM  REAL  ARRAY  DISCUSSED  ABOVE. 

C 

DIMENSION  B(2) 

COMMON  /CONI/  PI I,  P7,  P7TWO,  C22,  S22,  PI 2 

IW  IS  A  MACHINE  DEPENDENT  WRITE  DEVICE  NUMBER 

IW  =  I1MACH(2) 

IW=6 
C 

PII  =  4.*ATAN(1.) 

PIS  =  PII/8. 

P7  =  l./SQRT(2.) 

P7TWO  =  2.*P7 
C22  =  COS(PI8) 

S22  =  SIN(PI8) 

PI2  =  2.*PII 
N  =  1 

DO  10  1=1,15 
M  =  I 
N  =  N*2 

IF  (N.EQ.NFFT)  GO  TO  20 
10  CONTINUE 

WRITE  (IWf9999) 

9999  FORMAT  (30H  NFFT  NOT  A  POWER  OF  2  FOR  FFS) 

STOP 

20  CONTINUE 

B(2)  =  B(NFPT+1) 

DO  30  1=1, NFFT 

B(I)  =  B( I) /FLOAT (NFFT) 

30  CONTINUE 

DO  40  1=4, NFFT, 2 
B(I)  =  -B(I) 
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40  CONTINUE 

N8P0W  =  M/3 

REORDER  THE  INPUT  FOURIER  COEFFICIENTS 

CALL  0RD2(M,  B) 

CALL  ORDKM,  B) 

IF  (N8POW.EQ.0)  GO  TO  60 

PERFORM  THE  RADIX  8  ITERATIONS 

NN  =  N 

DO  50  IT=l,N8POW 
INT  =  N/NN 

CALL  R8SYN(INT,  NN,  B,  B(INT+1) ,  B(2*INT+1),  B(3*INT+1), 

*  B(4*INT+1),  B(5*INT+1),  B(6*INT+1) ,  B(7*INT+1) ,  B(l) , 

*  B(INT+1),  B(2*INT+1),  B(3*INT+1),  B(4*INT+1), 
B(5*INT+1) , 

*  B(6*INT+1),  B(7*INT+1)) 

NN  =  NN/8 

50  CONTINUE 

DO  A  RADIX  2  OR  RADIX  4  ITERATION  IF  ONE  IS  REQUIRED 

60  IF  (M-N8POW*3-l)  90,  80,  70 
70  INT  =  N/4 

CALL  R4SyN(INT,  B(l) ,  B(INT+1) ,  B(2*INT+1) ,  B(3*INT+1)) 

GO  TO  90 
80  INT  =  N/2 

CALL  R2TR(INT,  B(l) ,  B(INT+1)) 

90  RETURN 
END 


SUBROUTINE;  R2TR 

RADIX  2  ITERATION  SUBROUTINE 


SUBROUTINE  R2TR(INT,  B0,  Bl) 
DIMENSION  B0(2) ,  Bl(2) 

DO  10  K=1,INT 

T  =  B0(K)  +  B1(K) 

B1(K)  =  B0(K)  -  B1(K) 

B0(K)  =  T 
10  CONTINUE 
RETURN 
END 


C  SUBROUTINE;  R4TR 
C  RADIX  4  ITERATION  SUBROUTINE 


C 

SUBROUTINE  R4TR(INT,  B0 ,  Bl,  B2,  B3) 

DIMENSION  B0(2),  Bl(2),  B2(2) ,  B3(2) 

DO  10  K=1,INT 

R0  =  B0(K)  +  B2(K) 

R1  =  B1(K)  +  B3(K) 

B2(K)  s  B0(K)  -  B2(R) 

B3(K)  *  B1(K)  -  B3(K) 

B0(K)  *  R0  +  R1 
B1(K)  =  R0  -  R1 
10  CONTINUE 
RETURN 
END 
C 

C - 

C  SUBROUTINE;  R8TR 
C  RADIX  8  ITERATION  SUBROUTINE 

C - 

C 

SUBROUTINE  R8TR(INT,  NN,  BR0,  BRl ,  BR2,  BR3,  BR4,  BR5,  BR6 , 

BR7, 

*  BI0,  BIl,  BI2,  BI3,  BI4,  BIS,  BI6,  BI7> 

DIMENSION  L(15),  BR0(2),  BR1(2) ,  BR2(2),  BR3(2),  BR4(2), 

BR5(2)  , 

*  BR6(2)r  BR7(2),  BI0(2)  ,  BI1(2),  B12(2)  ,  BI3(2)  ,  BI4(2)  , 

*  BI5(2)  ,  BI6(2)  ,  BI7(2) 

COMMON  /CON/  PII,  P7,  P7TWO,  C22,  S22,  PI2 

EQUIVALENCE  (L15,L(1)),  (L14,L{2)),  (L13,L(3)),  (L12,L(4)), 

*  (L11,L(5)),  (L10,L(6)),  (L9,L(7)),  (L8,L(8)),  (L7,L(9)) 

*  (L6,L(10)),  (L5,L(11)),  (L4,L(12)),  (L3,L(13)), 

(L2,L(14)) , 

*  (L1,L(15)) 

C 

C  SET  OP  COUNTERS  SUCH  THAT  JTHET  STEPS  THROUGH  THE  ARGUMENTS 
C  OF  W,  JR  STEPS  THROUGH  STARTING  LOCATIONS  FOR  THE  REAL  PART  OF 
THE 

C  INTERMEDIATE  RESULTS  AND  JI  STEPS  THROUGH  STARTING  LOCATIONS 
C  OF  THE  IMAGINARY  PART  OP  THE  INTERMEDIATE  RESULTS. 

C 

L(l)  =  NN/8 
DO  40  K=2,15 

IP  (L(K-l)-2)  10,  20,  30 
10  L(K-l)  =  2 

20  L(K)  =  2 

GO  TO  40 

30  L(K)  =  L(K-l)/2 

40  CONTINUE 

PIOVN  *  PII/PLOAT(NN) 


w 


JI 

-  3 

JL 

»  2 

JR 

=  2 

DO 

120 

J1=2,L1, 

2 

DO 

120 

J2=J1,L2 

fLl 

DO 

120 

J3=J2,L3 

,L2 

DO 

120 

J4sJ3,L4 

rL3 

DO 

120 

J5=J4,L5 

rL4 

DO 

120 

J6»J5,L6 

,L5 

DO 

120 

J7=J6,L7 

,L6 

DO 

120 

J8=J7,L8 

DO 

120 

J9=J8,L9 

,L8 

DO 

120 

J10=J9,L10,: 

L9 

DO 

120 

J11*J10, 

Lll 

,L10 

DO 

120 

J12»J11, 

LI  2 

,L11 

DO 

120 

J13=J12, 

L13 

,L12 

DO 

120 

J14*J13, 

L14 

,L13 

DO 

120 

JTHET=J14,L15,L14 

TH2  =  JTHET  -  2 
IF  (TH2)  50,  50,  90 
50  DO  60  K=1,INT 

T0  =  BR0(K)  +  BR4(K) 

T1  =  BRl(K)  +  BR5(K) 

T2  =  BR2(K)  +  BR6(K) 

T3  =  BR3(K)  +  BR7(K) 

T4  =  BR0(K)  -  BR4(K) 

T5  =  BRl(K)  -  BR5(K) 

T6  =  BR2(K)  -  BR6(K) 

T7  =  BR3(K)  -  BR7(K) 
BR2(K)  =  T0  -  T2 

BR3(K)  =  T1  -  T3 

T0  =  T0  +  T2 
T1  =  Tl  +  T3 
BR0(K)  =  T0  +  Tl 

BRl(K)  =  T0  -  Tl 

PR  =  P7*{T5-T7) 

PI  =  P7*(T5+T7) 

BR4(K)  *  T4  +  PR 

BR7(K)  =  T6  +  PI 

BR6(K)  =  T4  -  PR 

BR5(K)  *  PI  -  T6 

60  CONTINUE 

IF  (NN-8)  120,  120,  70 
70  K0  =  INT*8  +  1 

KL  =  K0  +  INT  -  1 
DO  80  K»K0,KL 

PR  »  P7*(BI2(K)-BI6(K) ) 
PI  =  P7*(BI2(K)+BI6(K)) 
TR0  =  BI0(K)  +  PR 

TI0  -  BI4(K)  +  PI 

TR2  »  BI0(K)  -  PR 

TI2  »  BI4(K)  -  PI 

PR  =  P7*(BI3(K)-BI7(K)) 
PI  -  P7*(BI3(K)+BI7(K)) 
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TRl 

=  BIl(K)  +  PR 

Til 

=  BI5(K)  +  PI 

TR3 

»  BIl(K)  -  PR 

TI3 

=  BI5(K)  -  PI 

PR 

*  TR1*C22  -  TI1*S22 

PI 

=  TI1*C22  +  TR1*S22 

BI0(K)  =  TR0  +  PR 
BI6(K)  =  TR0  -  PR 
BI7(K)  =  TI0  +  PI 
BIl(K)  =  PI  -  TI0 

PR 

=  -TR3*S22  -  TI3*C22 

PI 

=  TR3*C22  -  TI3*S22 

BI2(K)  =  TR2  +  PR 
BI4(K)  =  TR2  -  PR 
BI5(K)  =  TI2  +  PI 
BI3(K)  =  PI  -  TI2 

CONTINUE 

GO  TO 

120 

ARG  = 

TH2*PIOVN 

Cl  = 

COS (ARG) 

SI  = 

SIN(ARG) 

C2  = 

Cl**2  -  Sl**2 

S2  = 

C1*S1  +  C1*S1 

C3  = 

C1*C2  -  S1*S2 

S3  = 

C2*S1  +  S2*C1 

C4  = 

C2**2  -  S2**2 

S4  = 

C2*S2  +  C2*S2 

C5  * 

C2*C3  -  S2*S3 

S5  = 

C3*S2  +  S3*C2 

C6  = 

C3**2  -  S3**2 

S6  = 

C3*S3  +  C3*S3 

C7  = 

C3*C4  -  S3*S4 

SI  = 

C4*S3  +  S4*C3 

INT8 

=  INT*8 

J0  ® 

JR*INT8  +  1 

K0  = 

JI*INT8  +  1 

JLAST 

=  J0  +  INT  -  ] 

DO  100  J=J0, JLAST 

K  = 

K0  +  J  -  J0 

TRl 

=  BR1(J)*C1  - 

BI1(K)*S1 

Til 

=  BR1(J)*S1  + 

BI1(K)*C1 

TR2 

*  BR2(J)*C2  - 

BI2(K) *S2 

TI2 

=  BR2(J)*S2  + 

BI2(K)*C2 

TR3 

=  BR3(J)*C3  - 

BI3(K)*S3 

TI3 

=  BR3(J)*S3  + 

BI3(K) *C3 

TR4 

=  BR4(J)*C4  - 

BI4(K) *S4 

TI4 

=  BR4(J)*S4  + 

BI4(K) *C4 

TR5 

=  BR5(J)*C5  - 

BI5(K) *S5 

TI5 

=  BR5(J)*S5  + 

BI5(K)*C5 

TR6 

=  BR6(J)*C6  - 

BI6(K) *S6 

TI6 

=  BR6{J)*S6  + 

BI6(K)*C6 

TR7 

«  BR7(J)*C7  - 

BI7(K) *S7 

TI7 

»  BR7(J) *S7  + 

BI7(K) *C7 

%  .  % 


t 


T0  =  BR0(J)  +  TR4 
T1  »  BI0(K)  +  TI4 

ID  TR4  »  BR0(J)  -  TR4 

f  TI4  *  BI0(K)  -  TI4 

T2  =  TRl  +  TR5 

T3  =  Til  +  TI5 

TR5  =  TRl  -  TR5 

TI5  »  Til  “  TI5 

M  T4  =  TR2  +  TR6 

T5  »  TI2  +  TI6 
TR6  =  TR2  -  TR6 

TI6  =  TI2  -  TI6 

T6  =  TR3  +  TR7 

T7  =  TI3  +  TI7 

TR7  *  TR3  -  TR7 

TI7  =  TI3  -  TI7 


C 


TR0  =  T0 

+  T4 

TI0  =  T1 

+  T5 

TR2  =  T0 

-  T4 

r 

TI2  *  T1 

-  T5 

TRl  «  T2 

+  T6 

Til  =  T3 

+  T7 

TR3  =  T2 

-  T6 

TI3  =  T3 

-  T7 

T0  =  TR4 

-  TI6 

k 

T1  *  TI4 

+  TR6 

T4  =  TR4 

+  TI6 

T5  *  TI4 

-  TR6 

T2  *  TR5 

“  TI7 

T3  =  TI5 

+  TR7 

T6  =  TR5 

+  TI7 

fi 

T7  =  TI5 

-  TR7 

BR0(J)  s 

TR0  +  TRl 

BI7(K)  * 

TI0  +  Til 

BI6(K)  » 

TR0  -  TRl 

BRl(J)  = 

Til  -  TI0 

' 

BR2(J)  = 

TR2  -  TI3 

BI5(K)  = 

TI2  +  TR3 

BI4(K}  s 

TR2  +  TI3 

V 

BR3(J)  = 

TR3  -  TI2 

PR  =  P7* 

(T2-T3) 

PI  =  P7* 

(T2+T3) 

BR4(J)  « 

T0  +  PR 

BI3(K)  * 

T1  +  PI 

BI2(K)  » 

T0  -  PR 

V' 

BR5(J)  = 

PI  -  T1 

•  • 

PR  =  -P7 

*(T6+T7) 

PI  *  P7* 

(T6-T7) . 

BR6(J)  > 

T4  +  PR 

L 

BIl(K)  * 

T5  +  PI 

BI0(K)  > 

T4  -  PR 

a' 

BR7(J)  - 

PI  -  T5 

L 
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SUBROUTINE  R4SYN(INT,  B0,  Bl,  B2,  B3) 
DIMENSION  B0(2),  Bl(2),  B2(2)r  B3(2} 
DO  10  K=1,INT 

T0  =  B0(K)  +  B1(K) 

T1  =  B0(K)  -  B1(K) 

T2  =  B2(K)  +  B2(K) 

T3  =  B3(K)  +  B3(K) 

B0(K)  =  T0  +  T2 

B2{K)  =  T0  -  T2 

B1(K)  =  T1  +  T3 

B3(K)  *  T1  -  T3 

CONTINUE 
RETURN 


C  SUBROUTINE:  R8SYN 
C  RADIX  8  SYNTHESIS  SUBROUTINE 


SUBROUTINE  R8SYN(INT,  NN,  BR0,  BRl ,  BR2,  BR3 ,  BR4,  BR5,  BR6 

BR7 

BI0,  BIl,  BI2,  BI3,  BI4,  BIS,  BI6,  BI7) 

DIMENSION  L(15),  BR0(2),  BRl(2),  BR2(2),  BR3(2),  BR4(2), 
BR5(2) , 

*  BR6(2),  BR7(2),  BI0(2),  BI1(2),  BI2(2),  BI3(2),  BI4(2), 

*  BI5(2) ,  BI6(2) ,  BI7(2) 

COMMON  /CONI/  PI I,  P7,  P7TWO,  C22,  S22,  PI 2 

EQUIVALENCE  (L15,L(1)),  (L14,L(2)),  (L13,L(3)),  (L12,L(4)), 

*  (L11,L{5)),  (L10,L(6)),  (L9,L(7)),  (L8,L(8)),  (L7,L(9)) 

*  (L6,L(10)),  (L5,L(11)),  (L4,L(12)),  (L3,L(13)), 

(L2,L(14)) , 

*  (L1,L(15)) 

L(l)  =  NN/8 
DO  40  K=2,15 


IF  (L(K-l)-2)  10,  20,  30 
10  L(K-l)  »  2 

20  L(K)  =  2 

GO  TO  40 

30  L(K)  =  L(K-l)/2 

40  CONTINUE 

PIOVN  =  PII/PLOAT(NN) 

JI  =  3 
JL  =  2 
JR  =  2 


DO 

120 

J1=2,L1, 

2 

DO 

120 

J2=J1,L2 

/LI 

DO 

120 

J3=J2,L3 

,L2 

DO 

120 

J4=J3,L4 

/L3 

DO 

120 

J5=J4,L5 

/L4 

DO 

120 

J6=J5,L6 

/L5 

DO 

120 

J7=J6,L7 

,L6 

DO 

120 

J8=J7,L8 

/L7 

DO 

120 

J9=J8,L9 

,L8 

DO 

120 

J10=J9,L10,: 

L9 

DO 

120 

J11=J10, 

Lll 

,L10 

DO 

120 

J12=J11, 

LI  2 

,L11 

DO 

120 

J13=J12, 

L13 

,L12 

DO 

120 

J14=J13, 

LI  4 

,L13 

DO 

120 

JTHET=J14,L15,L14 

TH2  =  JTHET  -  2 
IP  (TH2)  50,  50,  90 
50  DO  60  K=1,INT 

T0  =  BR0(K)  +  BRl(K) 
Tl  =  BR0(K)  -  BRl(K) 
T2  =  BR2(K)  +  BR2(K) 
T3  =  BR3(K)  +  BR3(K) 
T4  =  BR4(K)  +  BR6(K) 
T6  =  BR7(K)  -  BR5(K) 
T5  =  BR4(K)  -  BR6(K) 
T7  =  BR7(K)  +  BR5(K) 
PR  =  P7*(T7+T5) 

PI  *  P7*(T7-T5) 

TT0  =  T0  +  T2 
TTl  *  Tl  +  T3 
T2  =  T0  -  T2 
T3  *  Tl  -  T3 
T4  “  T4  +  T4 
T5  =  PR  +  PR 
T6  =  T6  +  T6 
T7  »  PI  +  PI 
BR0(K)  =  TT0  +  T4 
BRl(K)  =  TTl  +  T5 
BR2(K)  *  T2  +  T6 
BR3(K)  =  T3  +  T7 
BR4(K)  «  TT0  -  T4 
BR5(K)  =  TTl  -  T5 
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BR7(K)  =  T3  -  T7 
CONTINUE 

IF  (NN-8)  120,  120,  70 
K0  «  INT*8  +  1 
KL  »  K0  +  INT  -  1 
DO  80  K«K0,KL 

T1  >  BI0(K)  +  BI6(K) 

T2  «  BI7(K)  -  BIl(K) 

T3  »  BI0(K)  -  BI6(K) 

T4  »  BI7(K)  +  BIl(K) 

PR  »  T3*C22  +  T4*S22 
PI  =  T4*C22  -  T3*S22 
T5  =  BI2(K)  +  BI4(K) 

T6  =  BI5(K)  -  BI3(K) 

T7  »  BI2(K)  -  BI4(K) 

T8  =  BI5(K)  +  BI3(K) 

RR  =  T8*C22  -  T7*S22 
RI  *  -T8*S22  -  T7*C22 
BI0(K)  *  (T1+T5)  +  (T1+T5) 
BI4(K)  =  (T2+T6)  +  (T2+T6) 
BIl(K)  =  (PR+RR)  +  (PR+RR) 
BI5(K)  =  (PI+RI)  +  (PI+RI) 
T5  =  T1  -  T5 
T6  =  T2  -  T6 
BI2(K)  =  P7TWO*(T6+T5) 
BI6(K)  »  P7TWO*(T6-T5) 

RR  *  PR  -  RR 
RI  =  PI  “  RI 
BI3(K)  =  P7TWO*(RI+RR) 
BI7(K)  =  P7TWO*(RI-RR) 
CONTINUE 
GO  TO  120 
ARG  »  TH2*PIOVN 


Cl 

= 

cos (ARG) 

SI 

s 

-SIN (ARG) 

C2 

s 

Cl**2 

- 

Sl**2 

S2 

s 

C1*S1 

+ 

C1*S1 

C3 

s 

C1*C2 

- 

S1*S2 

S3 

s 

C2*S1 

+ 

S2*C1 

C4 

s 

C2**2 

- 

S2**2 

S4 

St 

C2*S2 

+ 

C2*S2 

C5 

s 

C2*C3 

- 

S2*S3 

S5 

s 

C3*S2 

+ 

S3*C2 

C6 

s 

C3**2 

- 

S3**2 

S6 

s 

C3*S3 

+ 

C3*S3 

C7 

s 

C3*C4 

- 

S3*S4 

S7 

s 

C4*S3 

+ 

S4*C3 

INT8 

=  INT* 8 

J0 

ss 

JR*INT8 

+  1 

K0 

a 

JI*INT8 

+  1 

JLAST  »  J0 

+ 

INT  - 

DO  100  JsJ0,JLAST 


K  *  K0  +  J  -  J0 
TR0  =  BR0(J)  +  BI6(K) 
TI0  =  BI7(K)  -  BRl(J) 
TRl  =  BR0(J)  -  BI6(K) 
Til  =  BI7(K)  +  BRl(J) 
TR2  =  BR2(J)  +  BI4(K) 
TI2  =  BI5(K)  -  BR3(J) 
TR3  =  BI5(K)  +  BR3(J) 
TI3  =  BI4(K)  -  BR2(J) 
TR4  =  BR4(J)  +  BI2(K) 
TI4  »  BI3(K)  -  BR5(J) 
T0  =  BR4(J)  -  BI2(K) 
T1  =  BI3(K)  +  BR5(J) 
TR5  =  P7*(T0+T1) 

T15  =  P7*(T1-T0) 

TR6  =  BR6(J)  +  BI0(K) 
TI6  =  BIl(K)  -  BR7(J) 
T0  =  BR6(J)  -  BI0(K) 
T1  =  BIl(K)  +  BR7(J) 
TR7  =  “P7*(T0-T1) 

T17  *  -P7*(T1+T0) 

T0  =  TR0  +  TR2 
T1  =  TI0  +  T12 
T2  =  TRl  +  TR3 
T3  =  Til  +  TI3 
TR2  =  TR0  -  TR2 
TI2  =  TI0  -  TI2 
TR3  =  TRl  -  TR3 
TI3  »  Til  -  TI3 
T4  =  TR4  +  TR6 
T5  =  TI4  +  TI6 
T6  =  TR5  +  TR7 
T7  =  TI5  +  TI7 
TTR6  =  TI4  -  TI6 
TI6  =  TR6  -  TR4 
TTR7  =  TI5  -  TI7 
TI7  =  TR7  -  TR5 


BR0(J) 
BI0(K) 
BRl(J) 
BIl(K) 
BR2( J) 
BI2(K) 
BR3(J) 
BI3(K) 
BR4(J) 
BI4(K) 
BR5(J) 
BI5(K) 
BR6(J) 
BI6(K) 
BR7(J) 
BI7(K) 


T0  +  T4 
T1  +  T5 
C1*(T2+T6)  - 
C1*(T3+T7)  + 
C2*(TR2+TTR6) 
C2*(TI2+TI6) 
C3*(TR3+TTR7) 
C3*(TI3+TI7) 
C4*(T0-T4)  - 
C4*(T1-T5)  + 
C5*(T2-T6)  - 
C5*(T3-T7)  + 
C6*(TR2-TTR6) 
C6*(TI2-TI6) 
C7*(TR3-TTR7) 
C7*(TI3-TI7) 


S1*(T3+T7) 

S1*(T2+T6) 

-  S2*(TI2+TI6) 
+  S2*(TR2+TTR6) 

-  S3*(TI3+TI7) 
+  S3*(TR3+TTR7) 
S4*(T1-T5) 
S4*(T0-T4) 
S5*(T3-T7) 
S5*(T2-T6) 

-  S6*(TI2-TI6) 
+  S6*(TR2-TTR6) 

-  S7*(TI3-TI7) 
+  S7*(TR3-TTR7) 


100  CONTINUE 

JR  =  JR  +  2 

JI  =  JI  -  2 

IP  (JI-JL)  110,  110,  120 
110  JI  »  2*JR  -  1 

JL  »  JR 
120  CONTINUE 
RETURN 
END 
C 


C  SUBROUTINE :  ORDl 
C  IN-PLACE  REORDERING  SUBROUTINE 


C  SUBROUTINE :  ORD2 
C  IN-PLACE  REORDERING  SUBROUTINE 

C - 

C 

SUBROUTINE  ORD2(M,  B) 

DIMENSION  L(15) ,  B(2) 

EQUIVALENCE  (L15,L(1)),  (L14,L(2)),  (L13,L(3)),  (L12,L(4)), 

*  (L11,L(5)),  (L10,L(6)),  (L9,L(7)),  (L8,L(8)),  (L7,L(9)) 

*  (L6,L(10)),  (L5,L(11)),  (L4,L(12)),  (L3,L(13)), 

(L2,L(14)) , 

*  (L1,L(15)) 

N  -  2**M 

L(l)  -  N 
DO  10  K»2,M 

L(K)  »  L(K-l)/2 


10 

CONTINUE 

DO 

20 

K«M,14 

L(K+1)  *  2 

20 

CONTINUE 

IJ 

■  : 

2 

DO 

40 

J1*2,L1,2 

DO 

40 

J2»J1,L2,L1 

DO 

40 

J3«J2,L3,L2 

DO 

40 

J4>J3,L4,L3 

DO 

40 

J5sJ4,L5,L4 

DO 

40 

J6»J5,L6,L5 

DO 

40 

J7-J6,L7,L6 

DO 

40 

J8»J7,L8,L7 

DO 

40 

J9»J8,L9,L8 

DO 

40 

J10aJ9rL10,L9 

DO 

40 

J11=J10,L11,L10 

DO 

40 

J12=J11,L12,L11 

DO 

40 

J13»J12,L13,L12 

DO 

40 

J14«J13,L14,L13 

DO 

40 

JI=J14,L15,L14 

IP 

(IJ-JI)  30,  40,  40 

30 

1 

r  = 

B(IJ-l) 

B(IJ-l)  =  B(JI-l) 
B(JI-l)  *  T 

r  = 

B(IJ) 

B(IJ)  =  B(JI) 

B(JI)  =  T 

40 

IJ  = 

=  IJ  +  2 

RETURN 

C  SUBROUTINE:  FFT842 
C  FAST  FOURIER  TRANSFORM  FOR  N=2**M 
C  COMPLEX  INPUT 


SUBROUTINE  PFT842(IN,  N,  X,  Y) 


C  THIS  PROGRAM  REPLACES  THE  VECTOR  Z=X+iy  BY  ITS  FINITE 
C  DISCRETE,  COMPLEX  FOURIER  TRANSFORM  IF  IN=0.  THE  INVERSE 
TRANSFORM 

C  IS  CALCULATED  FOR  IN=1.  IT  PERFORMS  AS  MANY  BASE 
C  8  ITERATIONS  AS  POSSIBLE  AND  THEN  FINISHES  WITH  A  BASE  4 
ITERATION 

C  OR  A  BASE  2  ITERATION  IF  NEEDED. 


C  THE  SUBROUTINE  IS  CALLED  AS  SUBROUTINE  FFT842  (IN,N,X,Y). 

C  THE  INTEGER  N  (A  POWER  OF  2) ,  THE  N  REAL  LOCATION  ARRAY  X,  AND 
C  THE  N  REAL  LOCATION  ARRAY  Y  MUST  BE  SUPPLIED  TO  THE  SUBROUTINE. 


on  non  non  \o  o  noono 


DIMENSION  X(2),  Y(2),  L(15) 

COMMON  /C0N2/  PI 2,  P7 

EQUIVALENCE  (L15,L(1)),  (L14,L(2)),  {L13,L(3)),  (L12,L(4)), 

*  (L11,L(5)),  (L10,L(6)),  (L9,L(7)),  (L8,L(8)),  {L7,L(9)), 

*  (L6,L(10)),  (L5,L(11)),  (L4,L(12)),  (L3,L(13)), 

(L2,L(14)) , 

*  (L1,L(15)) 


IW  IS  A  MACHINE  DEPENDENT  WRITE  DEVICE  NUMBER 

IW  =  I1MACH(2) 

IW=6 

PI2  =  8.*ATAN(1.) 

P7  =  l./SQRT(2.) 

DO  10  1=1,15 
M  =  I 
NT  =  2**1 

IF  (N.EQ.NT)  GO  TO  20 
10  CONTINUE 

WRITE  (IW,9999) 

199  FORMAT  (35H  N  IS  NOT  A  POWER  OF  TWO  FOR  FFT842) 

STOP 

20  N2P0W  =  M 
NTHPO  =  N 
FN  =  NTHPO 

IF  (IN.EQ.l)  GO  TO  40 
DO  30  1=1, NTHPO 
Y(I)  =  -Y(I) 

30  CONTINUE 
40  N8POW  =  N2POW/3 

IF  (N8POW.EQ.0)  GO  TO  60 

RADIX  8  PASSES, IF  ANY, 

DO  50  IPASS=1,N8P0W 

NXTLT  =  2**(N2POW-3*IPASS) 

LENGT  =  8*NXTLT 

CALL  R8TX(NXTLT,  NTHPO,  LENGT,  X(l),  X(NXTLT+1), 

;(2*NXTLT+1)  , 

*  X(3*NXTLT+1) ,  X{4*NXTLT+1) ,  X ( 5*NXTLT+1) ,  X ( 6*NXTLT+1) , 

*  X(7*NXTLT+1) ,  Y(l),  Y(NXTLT+1),  Y ( 2*NXTLT+1) , 
(3*NXTLT+1) , 

*  Y{4*NXTLT+1) ,  Y(5*NXTLT+1) ,  Y ( 6*NXTLT+1) ,  Y ( 7*NXTLT+1) ) 
50  CONTINUE 

IS  THERE  A  FOUR  FACTOR  LEFT 

60  IF  (N2POW-3*N8POW-l)  90,  70,  80 

GO  THROUGH  THE  BASE  2  ITERATION 


B-38 


non 


70  CALL  R2TX(NTHP0,  X(l) ,  X(2),  Y(l) ,  y(2)) 

GO  TO  90 

GO  THROUGH  THE  BASE  4  ITERATION 

80  CALL  R4TX(NTHPOr  X(l),  X(2)r  X(3)f  X(4),  Yd),  Y(2)  ,  Y(3) 
Y(4)) 

90  DO  110  J=l,15 
L(J)  =  1 

IF  (J-N2POW)  100,  100,  110 
100  L(J)  =  2**(N2POW+l-J) 

110  CONTINUE 
IJ  =  1 

DO  130  J1=1,L1 
DO  130  J2=J1,L2,L1 
DO  130  J3=J2,L3,L2 
DO  130  J4=J3,L4,L3 
DO  130  J5=J4,L5,L4 
DO  130  J6=J5,L6,L5 
DO  130  J7=J6,L7,L6 
DO  130  J8=J7,L8,L7 
DO  130  J9=J8,L9,L8 
DO  130  J10=J9,L10,L9 
DO  130  J11=J10,L11,L10 
DO  130  J12=J11,L12,L11 
DO  130  J13=J12,L13,L12 
DO  130  J14=J13,L14,L13 
DO  130  JI=J14,L15,L14 
IF  (IJ-JI)  120,  130,  130 
120  R  =  X(IJ) 

X(IJ)  =  X(JI) 

X(JI)  =  R 
FI  =  Y(IJ) 

Y(IJ)  =  Y(JI) 

Y(JI)  =  FI 
130  IJ  =  IJ  +  1 

IF  (IN.EQ.l)  GO  TO  150 
DO  140  I=l,NTHPO 
Y(I)  =  -Y(I) 

140  CONTINUE 
GO  TO  170 

150  DO  160  I=l,NTHPO 
X(I)  =  X(I)/FN 
Y(I)  =  Y(I)/FN 
160  CONTINUE 
170  RETURN 


C  SUBROUTINE:  R2TX 
C  RADIX  2  ITERATION  SUBROUTINE 

C - 


C 

SUBROUTINE  R2TX(NTHPO,  CR0,  CRl,  CI0,  CIl) 

DIMENSION  CR0(2),  CR1(2),  CI0(2) ,  CI1(2) 

DO  10  K«l,NTHPO,2 
R1  *  CR0(K)  +  CRl(K) 

CR1{K)  »  CR0{K)  -  CRl(K) 

CR0(K)  =  R1 
FIl  =  CI0(K)  +  CIl(K) 

CIl(K)  =  CI0(K)  -  CIl(K) 

CI0(K)  =  FIl 
10  CONTINUE 
RETURN 
END 

C 

C - 

C  SUBROUTINE :  R4TX 

C  RADIX  4  ITERATION  SUBROUTINE 

C - 

C 

SUBROUTINE  R4TX(NTHPO,  CR0 ,  CRl,  CR2,  CR3,  CI0,  CIl,  CI2 

CI3) 

DIMENSION  CR0(2),  CR1(2),  CR2(2),  CR3(2),  CI0(2) ,  CIl (2) 
CI2(2)  , 

*  CI3(2) 

DO  10  K=l,NTHPO,4 
R1  =  CR0(K)  +  CR2(K) 

R2  =  CR0(K)  -  CR2(K) 

R3  *  CRl(K)  +  CR3(K} 

R4  =  CRl(K)  -  CR3(K) 

FIl  =  CI0(K)  +  CI2(K) 

FI2  =  CI0(K)  -  CI2(K) 

FI3  =  CIl(K)  +  CI3(K) 

FI4  =  CIl(K)  -  CI3(K) 

CR0(K)  =  R1  +  R3 
CI0(K)  =  FIl  +  FI3 
CRl(K)  =  R1  -  R3 
CIl(K)  =  FIl  -  FI3 
CR2(K)  =  R2  -  FI4 
CI2(K)  =  FI2  +  R4 
CR3(K)  =  R2  +  FI4 
CI3(K)  *  FI2  -  R4 
10  CONTINUE 
RETURN 
END 


C  SUBROUTINE;  R8TX 
C  RADIX  8  ITERATION  SUBROUTINE 

C - 

c 

SUBROUTINE  R8TX(NXTLT,  NTHPO,  LENGT,  CR0,  CRl ,  CR2,  CR3,  CR4, 

*  CR5,  CR6,  CR7,  CI0,  CIl,  CI2,  CI3,  CI4,  CIS,  CI6,  CI7) 
DIMENSION  CR0(2),  CRl(2),  CR2(2) ,  CR3(2),  CR4(2),  CR5(2), 

CR6(2)  , 

*  CR7(2),  CI1(2),  CI2(2),  CI3{2)  ,  CI4(2),  CI5(2),  CI6(2), 

*  CI7(2) ,  CI0(2) 

COMMON  /CON2/  PI 2,  P7 

C 

SCALE  =  PI 2/PLOAT( LENGT) 

DO  30  J*1,NXTLT 

ARG  =  FLOAT { J-1) * SCALE 


Cl 

= 

cos (ARG) 

SI 

=: 

SIN (ARG) 

C2 

= 

Cl**2 

- 

Sl**2 

S2 

s 

C1*S1 

+ 

C1*S1 

C3 

= 

C1*C2 

- 

S1*S2 

S3 

= 

C2*S1 

+ 

S2*C1 

C4 

= 

C2**2 

- 

S2**2 

S4 

s: 

C2*S2 

+ 

C2*S2 

C5 

= 

C2*C3 

- 

S2*S3 

S5 

=: 

C3*S2 

+ 

S3*C2 

C6 

s 

C3**2 

- 

S3**2 

S6 

s 

C3*S3 

+ 

C3*S3 

C7 

s 

C3*C4 

- 

S3*S4 

S7 

s 

C4*S3 

+ 

S4*C3 

DO  20 

K= 

=J,NTHPO, LENGT 

AR0 

a 

CR0(K) 

+ 

CR4(K) 

ARl 

s 

CRl(K) 

+ 

CR5(K) 

AR2 

= 

CR2(K) 

+ 

CR6(K) 

AR3 

= 

CR3(K) 

+ 

CR7(K) 

AR4 

3= 

CR0(K) 

- 

CR4(K) 

AR5 

s 

CRl(K) 

- 

CR5(K) 

AR6 

= 

CR2(K) 

- 

CR6(K) 

AR7 

S 

CR3(K) 

- 

CR7(K) 

AI0 

3= 

CI0(K) 

+ 

CI4(K) 

All 

S 

CIl(K) 

+ 

CI5(K) 

AI2 

S 

CI2(K) 

+ 

CI6(K) 

AI3 

S 

CI3(K) 

+ 

CI7(K) 

AI4 

S 

CI0(K) 

- 

CI4(K) 

AI5 

as 

CIl(K) 

- 

CI5(K) 

AI6 

= 

CI2(K) 

- 

CI6(K) 

AI7 

s 

CI3(K) 

- 

CI7(K) 

BR0 

s 

AR0  + 

AR2 

BRl 

= 

ARl  + 

AR3 

BR2 

= 

AR0  - 

AR2 

BR3 

as 

ARl  - 

AR3 

BR4 

s 

AR4  - 

AI6 

BR5 

s 

AR5  - 

AI7 

BR6 

s 

AR4  + 

AI6 

-f 
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BR7  =  AR5  + 
BI0  =  AI0  + 
BIl  =  All  + 
BI2  =  AI0  - 
BI3  »  All  - 
BI4  =  AI4  + 
BIS  =  AI5  + 
BI6  =  AI4  - 
BI7  =  AI5  - 
CR0(K)  s  BR0 


AI7 
AI2 
AI3 
AI2 
AI3 
AR6 
AR7 
AR6 
AR7 
+  BRl 


CI0(K)  =  BI0  +  BIl 
IF  (J.LE.l)  GO  TO  10 
CRl(K)  =  C4*(BR0-BR1) 
CIl(K)  =  C4*(BI0-BI1) 
CR2{K)  =  C2*(BR2-BI3) 
CI2(K)  =  C2*(BI2+BR3) 
CR3(K)  =  C6*(BR2+BI3) 
CI3(K)  =  C6*(BI2-BR3) 
TR  =  P7*(BR5-BI5) 

TI  =  P7*(BR5+BI5) 
CR4(K)  =  C1*(BR4+TR) 
CI4(K)  =  C1*(BI4+TI) 
CR5(K)  =  C5*(BR4-TR) 
CI5(K)  =  C5*(BI4-TI) 
TR  =  >P7*(BR7+BI7) 

TI  =  P7*{BR7-BI7) 
CR6(K)  «  C3*(BR6+TR) 
CI6(K)  =  C3*(BI6+TI) 
CR7(K)  =  C7*(BR6-TR) 
CI7(K)  =  C7*(BI6-TI) 
GO  TO  20 

CRl(K)  =  BR0  -  BRl 
CIl(K)  »  BI0  -  BIl 
CR2(K)  =  BR2  -  BI3 
CI2(K)  =  BI2  +  BR3 
CR3(K)  =  BR2  +  BI3 
CI3(K)  =  BI2  -  BR3 
TR  =  P7*(BR5-BI5) 

TI  =  P7*(BR5+BI5) 
CR4(K)  =  BR4  +  TR 
CI4(K)  =  BI4  +  TI 
CR5(K)  =  BR4  -  TR 
CI5(K)  =  BI4  -  TI 
TR  =  -P7*(BR7+BI7) 

TI  =  P7*(BR7-BI7) 


S4*(BI0-BI1) 

S4*(BR0-BR1) 

S2*(BI2+BR3) 

S2*(BR2-BI3) 

S6*(BI2-BR3) 

S6*(BR2+BI3) 


-  S1*(BI4+TI) 
+  S1*(BR4+TR) 

-  S5*(BI4-TI) 
+  S5*(BR4-TR) 


-  S3*(BI6+TI) 
+  S3*(BR6+TR) 

-  S7*(BI6-TI) 
+  S7*(BR6-TR) 


CR6(K) 

CI6(K) 

CR7(K) 

CI7(K) 

CONTINUE 

CONTINUE 

RETURN 


=  BR6  + 
=  BI6  + 
=  BR6  “ 
=  BI6  - 


***************************************************************** 

CCXTRA.POR 

***************************************************************** 


SUBROUTINE  RP(X,Y,R,TH,N,ID) 

C  REC  -  POLAR  CONVERSION 
C  ID  NOT  =1  =»>  R->P 

C  ID  =1  ==>  P->R 

REAL*8  X(16384) ,Y(16384) ,R(16384) rTH(16384} 
REAL*8  DSQRT,DATAN2,DCOS,DSIN 
IP(ID.EQ.l)GO  TO  10 
DO  5  1=1, N 

R(I)=DSQRT(X(I)*X(I)+Y(I)*Y(I) ) 

TH(I)=0.D0 

IP(Y(I) .EQ.0.D0.AND.X(I) .EQ.0.D0)GO  TO  5 
TH(I)=DATAN2(Y(I) ,X(I) ) 

5  CONTINUE 
RETURN 
10  CONTINUE 
DO  15  1=1, N 
X(I)=R(I)*DCOS(TH(I)  ) 

15  Y(I) =R(I) *DSIN(TH(I)  ) 

RETURN 

END 


SUBROUTINE  MRIDP ( LOG2N,X ,Y ,SIGN) 

PORTRAN  VERSION 
MIXED  RADIX  POURIER  TRANSPORM 
INTEGER  LOG2N 
REAL*8  X(16384) ,Y(16384) 

DIMENSION  S(13) ,U(13) 

INTEGER  J1,J2,J3,J4,JT,N,M4 

REAL*8  ARG,C1,C2,C3,S1,S2,S3,R1,R2,R3,R4,R5,R6,R7,R8,T 
REAL*8  DCOS,DSIN 


INTEGER  A,B,C,D,E,P,G,H,I,J,K,L,M, 

1  BS,CS,DS,ES,PS,GS,HS,IS,JS,KS,LS,MS, 

2  AL,BL,CL,DL,EL,PL,GL,HL,IL,JL,KL,ML,S,U 
EQUIVALENCE 

(BS,S(2)) ,(CS,S(3)) ,(DS,S(4)) ,(ES,S(5)) ,(PS,S(6))  , 

1 

(GS,S(7)) ,{HS,S(8)) ,(IS,S(9)) ,(JS,S(10)) , (KS,S(11) ) , (LS,S(12) ) , 
2 

(MS,S{13)) ,(AL,U(1)) ,(BL,U(2)) ,(CL,U(3)) ,(DL,U(4)) ,(EL,U(5))  , 

3 

(PL,U(6)) ,(GL,U(7)) ,(HL,U(8)) , (IL,U(9) ) , ( JL,U(10) ) ,(KL,U(11)) , 

4  (LL,U(12)) ,(ML,U(13)) 


c 


N=2**LOG2N 

IF  (SIGN)  8000,8000,8002 

8000  DO  8001  I»1,N 

8001  y(I)=-Y(I) 

8002  CONTINUE 

IF  (LOG2N-1)  500,500,501 
501  CONTINUE 

DO  400  K»2,LOG2N,2 
M=2**(LOG2N-K) 

M4=4*M 

DO  300  J=1,M 

ARG=6 . 2  83 1 5D0  *DBLE ( FLOAT ( J-1 ) /FLOAT ( M4 ) ) 
Cl«DCOS(ARG) 

SlsDSIN(ARG) 

C2*C1*C1-S1*S1 

S2=C1*S1+C1*S1 

C3=C2*C1-S2*S1 

S3=C2*S1+S2*C1 

DO  200  I=M4,N,H4 

J1=I+J-M4 

J2=J1+M 

J3=J2+M 

J4=J3+M 

R1=X(J1) +X(J3) 

R2=X(J1)-X(J3) 

R3«Y(J1) +Y(J3) 

R4«Y(J1)-Y(J3) 

R5=X(J2)+X(J4) 

R6=X(J2)-X(J4) 

R7=Y(J2)+Y(J4) 

R8=Y(J2)-Y(J4) 

X(J1)=R1+R5 
Y(J1) =R3+R7 

IF(ARG)  101,100,101 
101  CONTINUE 

X(J3) =(R2+R8) *C1+(R4-R6) *S1 
Y( J3) =(R4-R6) *C1-(R2+R8) *S1 
X(J2)=(R1-R5)*C2+(R3-R7)*S2 
Y( J2) =(R3-R7) *C2-(R1-R5) *S2 
X(J4)  =  (R2-R8) *C3+(R4+R6)  *S3 
Y(J4) =(R4+R6) *C3-(R2-R8)  *S3 
GO  TO  200 
100  CONTINUE 

X(J3)=R2+R8 
Y(J3)-R4-R6 
X(J2)»R1-R5 
Y(J2)»R3-R7 
X(J4)*R2-R8 
Y(J4)»R4+R6 
200  CONTINUE 

300  CONTINUE 

400  CONTINUE 


500  CONTINUE 

ITEST=L0G2N-(L0G2N/2*2) 
IF(ITEST)  701,700,701 
701  CONTINUE 

DO  600  I»1,N,2 
Rl*X(I)+X{I+l) 

R2«X(I)-X(I+1) 

R3=Y(I)+Y(I+1) 

R4»Y(I)-Y(I+1) 

X{I)»R1 
Y(I)=R3 
X(I+1) »R2 
Y(I+1) =R4 
600  CONTINUE 

700  CONTINUE 

MS*N/2 
NLsN 

DO  800  K=2,12 
J=14-K 
S(J)=1 
U(J)=S(J+1) 

IF(S(J+1)“1)  7701,7701,7700 

7700  S(J)=S(J+l)/2 

7701  CONTINUE 

800  CONTINUE 

AL*BS 

JT=0 

DO  900  A»1,AL 
DO  900  BsA,BL,BS 
DO  900  C»=B,CL,CS 
DO  900  D=C,DL,DS 
DO  900  E=D,EL,ES 
DO  900  F=E,FL,FS 
DO  900  G=F,GL,GS 
DO  900  H»G,HL,HS 
DO  900  I=H,IL,IS 
DO  900  J=I,JL,JS 
DO  900  K=J,KL,KS 
DO  900  L-K,LL,LS 
DO  900  H»L,ML,MS 
JT=JT+1 

IF(JT-M)  900,900,901 
901  CONTINUE 

T=X(JT) 

X(JT)«X(M) 

X(M) =T 
T=Y(JT) 

Y(JT)=Y(M) 

Y(M)»T 

900  CONTINUE 
RETURN 
END 


B-45 


SHIFTF.FOR 

***************************************************************** 


C  •  •  • 

C  •  •  • 
ZERO 
C  •  •  • 

C  •  •  • 

c. . . 

C  •  •  • 

c. . . 

* 

c. . . 

* 

c. .  . 

*  y 

c. . . 

* 

C  •  •  • 

* 

C  •  •  • 

* 

C  •  • « 

0 « « • 

C  •  •  • 

C  •  •  • 

C  •  •  • 

C  •  •  • 

C  •  •  • 

c. .  . 
c. .  . 

C  •  •  • 

c. . . 

0  •  •  • 

C  •  •  • 

C  •  •  • 

C  •  •  • 

c. . . 

0  •  •  • 

C  •  •  • 

C  •  •  • 

C  •  •  • 

C  •  •  • 


SUBROUTINE  SHIFTF(F,N) 

SHIFTS  PHASE  OF  ARRAY  SO  THAT  ANGLE  OF  FIRST  ELEMENT  IS 


GIVEN  THE  FIRST  ELEMENT  a+ib  AND  ANY  ELEMENT  X+iY 


sqrt(a**2+b**2)  * 


*  b  and 


angle 


*  ANGLE 


******************** 


******************** 


SHIFT  THE  PHASE  BY  angle 


*  Y1 

* 


ANGLE-angle  * 

************************** 


“  SIN{angle) ®b/sqrt (a**2+b**2) 

COS (angle) =a/sqrt (a**2+b**2) 

I  XI  I  I  COS(angle}  SIN(sngle) I  I  X  | 

I  I  »  I  III 

I  Y1  I  l-SIN(angle)  COS(angle) I  I  Y  1 

SUCH  THAT  t 

xi=  [a/sqrt(a**2+b**2) ]*X  +  [b/sqrt {a**2+b**2) ] *Y 
Yl— [b/sqrt(a**2+b**2)  ]  *X  +  [a/sqrt  (a**2+b**2)  ]  *Y 


(I******************************************************************* 

********* 

DIMENSION  F(l) 

C  •  •  • 

A-P(l) 

B«P(2) 


J  '•  *' •***•“  -  •  ■  * ‘  •*•>'**  •**^-*"  '•  .*•  •**  .*■  •  *-  .* ■  j*'  *'*  *  '  • 


Appendix  C 

PATTERN  RECOGNITION  SCATTER  PLOT  PRODUCTION 


Start  with  several  sampled  data  files  that  represent  the  data  from  which  pattern 
similarities  are  to  be  extracted.  Each  segment  of  data  to  be  analyzed  is  labeled.  Label 
files  must  be  created  with  the  ILS  command  $  LBF  0,  which  is  similar  to  the  manner 
in  which  sampled  data  files  are  created  with  the  $  FIL  0  command  as  in  this  example. 
Use  of  a  command  file  is  suggested  when  many  label  files  are  to  be  created. 

To  label  the  segments,  use  ILS  command  $  LBA.  Data  segments  must  have  been 
previously  selected  by  use  of  one  of  the  ILS  commands  $  DSP  or  $  CUR  so  that  this 
information  is  in  the  user’s  common  file.  A  label  file  contains  a  series  of  labels,  each 
one  of  which  is  a  two-line  record  of  ASCII  characters  that  describes  a  segment  of  in¬ 
terest  in  a  sampled  data  file.  Next,  use  ILS  command  $  QUR  to  extract  signal  features. 
Because  ILS  command  $  QUR  restricts  the  user  to  32  frequency  spectra,  a  local  version 
was  used  that  selects  32  evenly  spaced  spectra  from  the  output  of  a  2048-point  FFT. 
To  use  QUR  (or  XUR),  first  designate  the  label  file  that  contains  pointers  to  the  area  of 
interest  in  the  signal.  $  QUR  output  will  be  in  the  form  of  feature  record  files. 

Using  ILS  command  $  SME,  create  output  feature  record  files  which  contain  mean 
vector  and  eigen  records.  Since  $  SME  is  now  restricted  to  20  elements  in  analysis, 
first  extract  the  20  most  significant  elements  in  the  output  records  using  ILS  command 
$  MRE.  Use  the  same  elements  for  all  records  analyzed.  If  you  have  used  QUR  on  sin¬ 
gle  data  samples,  collect  the  samples  in  two  or  more  record  files  using  $  TRE  for  analy¬ 
sis  with  $  SME.  In  this  study,  ILS  command  $  SME  6  was  used.  One  output  record  file 
is  produced  by  $  SME  containing  the  mean  and  eigen  vector  records.  Next,  using  the 
files  input  to  $  SME  and  the  output  file  of  $  SME,  perform  a  principal  component  anal¬ 
ysis  of  those  input  files.  $  PCO  outputs  a  record  file  for  each  input  file.  Plot  the  results 
withS  PLR. 


Examples  are  given  on  the  following  pages. 


Example : 

Ten  sampled  data  files  exist  which  are  numbered  1405  through  1495 
in  steps  of  10. 

Create  the  label  file(s) . 

$  LBF  0  <cr> 

OPTIONS  TO  MANIPULATE  LABEL  FILES: 

ENTER  FILENAME  [ ,DK]  TO  SELECT  A  FILE 
ENTER  C  <RETURN>  TO  CREATE  A  FILE 
ENTER  I  <RETURN>  TO  SET  THE  INITIALS 

ENTER  <RETURN>  TO  EXIT 

->C 

ENTER  FILENAME  [ ,DK] 

->TAM1 

DRB0: [ILSMGR.DGNJTAMI.LAB  LABEL  DATA 

0  LABELS,  INITIALS: 

PRIMARY  LABEL  FILE 

OPTIONS  TO  MANIPULATE  LABEL  FILES: 

ENTER  FILENAME  [ ,DK]  TO  SELECT  A  FILE 
ENTER  C  <RETURN>  TO  CREATE  A  FILE 
ENTER  I  <RETURN>  TO  SET  THE  INITIALS 

ENTER  <RETURN>  TO  EXIT 

->I 

ENTER  THE  INITIALS 
->TAM 

DRB0: tILSMGR.DGN]TEST2.LAB  LABEL  DATA 

0  LABELS,  INITIALS:  TAM 
PRIMARY  LABEL  PILE 

OPTIONS  TO  MANIPULATE  LABEL  PILES: 

ENTER  FILENAME  [ ,DK]  TO  SELECT  A  PILE 
ENTER  C  <RETURN>  TO  CREATE  A  PILE 

ENTER  I  <RETURN>  TO  SET  THE  INITIALS 

ENTER  <RETURN>  TO  EXIT 

-><cr> 

In  like  manner  create  the  other  nine  label  files. 

Use  of  a  command  file  is  suggested  when  many  label  files  are  to 
be  created. 

Label  the  segments,  having  previously  indicated  the  range  by  the 
cursor  command  or  by  the  $  DSP  range. 

Point  to  the  label  file  by  issuing  the  command  $  LBF  TAMl . 

DRB0: IILSMGR.ILS]TAM1.LAB  LABEL  DATA 

1  LABELS,  INITIALS: 

PRIMARY  LABEL  FILE 
$  LBA 

FIELD-1  IS  OBTAINED  PROM  THE  FILE  HEADER,  @=ABORT 
ENTER  FIELD-2; FIELD-3; FIELD-4; FIELD-5; FIELD-6 
->TAMTEST 

List  the  label  file 
$  LLA 

DRB0: tlLSMGR.ILSJTAMl.LAB  22-MAY-84 


PAGE 


*****  label  1  ************************************************** 

TAM7EST  •  •  •  •  • 

1;  1;  1;2048;  2000 ;DRA0 : [ ILSMGR. ILS] WD1405 .  ; 27-JAN-84 ;TAM 

Label  the  rest  of  the  files. 

Extract  the  features 
$  XUR 

ENTER  MODE,  EXTRACTION  CODE  AND  NUMBER  OF  FEATURES 
MODE:  M  -  MEAN  FRAMES 

C  »  CONSECUTIVE  FRAMES 
FEATURE:  1  -  AUTOREGRESSIVE  COEFFICIENTS 

2  »  REFLECTION  COEFFICIENTS 

3  =  AUTOCORRELATION  COEFFICIENTS 

4  =  FREQUENCY  SPECTRA 

5  *  LP  CEPSTRAL  COEFFICIENTS 

6  =  MEL  CEPSTRAL  COEFFICIENTS 

7  »  RESONANCE  FREQ  AND  BAND  (PEAK  PICK) 

8  =  RESONANCE  FREQ  AND  BAND  (ROOT  SOLVE) 

9  =  NORMALIZED  RESONANCE  FREQ  AND  BAND  (ROOT  SOLVE) 

->4 

EXTRACTION  CODE  =  4 

USING  LABEL  FILE. ... DRB0 :[ ILSMGR. ILS] TXMl . LAB 
USING  RECORD  FILE. .. .DRB0 :[ ILSMGR. ILS] WD1000 . 

LEVEL  1 

FIELD-1  IS  OBTAINED  FROM  THE  FILE  HEADER,  @=ABORT 
ENTER  FIELD-2; FIELD-3; FIELD-4; FIELD-5; FIELD-6  FOR  SEARCH 
ENTER  <RETURN>  TO  START  SEARCH 
-> 

DRB0: (ILSMGR. ILS] WD1000.  RECORD  1  STORED 

And  so  on  with  the  rest  of  the  files. 

Using  $  MRE  move  the  first  twenty  elements  in  each  file 
to  another  file 
$  MRE 

MOVE  FEATURE  RECORDS  FROM  A  FILE  AND 

EXTRACT  SELECTED  ITEMS  AND  ELEMENTS  FROM  RECORDS 


ENTER. . .ELEMENT  NUMBER(S)  TO  EXTRACT  FROM  EACH  ITEM 
USE  <ALL>  <RETURN>  FOR  ALL  ELEMENTS 

USE  <F>  <RETURN>  TO  FINISH 

USE  <A>  < RETURN >  TO  ABORT  . 

MAXIMUM  OF  12  ENTRIES/LINE,  NEGATIVE  IMPLIES  INCLUSIVE 

->1,-20 

->P 

ENTER... ITEM  NUMBER(S)  TO  EXTRACT  FROM  EACH  RECORD 
USE  <ALL>  <RETURN>  FOR  ALL  ITEMS 
USE  <F>  <RETURN>  TO  FINISH 

USE  <A>  <RETURN>  TO  ABORT 

MAXIMUM  OF  12  ENTRIES/LINE,  NEGATIVE  IMPLIES  INCLUSIVE 

->ALL 

ENTER. . .RECORD  NUMBER(S)  OF  EACH  RECORD  TO  TRANSFER 
USE  <ALL>  <RETURN>  FOR  ALL  RECORDS 
USE  <F>  <RETURN>  TO  FINISH 

USE  <T>  <RETURN>  TO  TEST  ON  ELEMENTS  (MAX=10) 

USE  <SK>  MO.  PEC  <P.ETURn>  TO  SKIP  RECORDS 


USE  <AV>  NO.  REC  <RETURN>  TO  AVERAGE  RECORDS 

USE  <AVI>  NO,  REC  <RETURN>  TO  AVERAGE  ITEMS  IN  RECORDS 

USE  <A>  <RETURN>  TO  ABORT 

MAXIMUM  OF  12  ENTRIES/LINE,  NEGATIVE  IMPLIES  INCLUSIVE 

->ALL 

ENTER... PILE  NO.,  DISK  NO.  OF  INPUT  PILE 
->1010 

ENTER... FILE  NO.,  DISK  NO.  OF  OUTPUT  FILE 
->1021 


DRB0: [ILSMGR.ILS]WD1021. 

RECORD 

1 

STORED 

And  so  on  with  the  rest  of  the  files. 

Collect  the  records  in  two  files  for  analysis, 

» 

$  FIL  1021 

DRB0 : [ ILSMGR. ILS] WD1021 . 

RECORD 

DATA 

12  DK  BLKS,  1  RECORDS 

PRIMARY  FILE 
$  FIL  S1031 

DRB0: [ILSMGR. ILS] WD1031. 

DOES  NOT  EXIST 

SECONDARY  FILE 
$  OPN  S2 

$  TRE  1,1, 5, 1,2, 5 

USING  DRB0: [ILSMGR. ILS] WD1021. 
DRB0: [ILSMGR. ILS] WD1031. 

RECORD 

1 

STORED 

USING  DRB0: [ILSMGR. ILS] WD1022. 
DRB0: [ILSMGR. ILS] WD1031. 

RECORD 

2 

STORED 

USING  DRB0; [ILSMGR. ILS] WD1023. 
DRB0 : [ ILSMGR. ILS] WD1031 . 

RECORD 

3 

STORED 

USING  DRB0: [ILSMGR. ILS] WO1024. 
DRB0! [ILSMGR. ILS] WD1031. 

RECORD 

4 

STORED 

USING  DRB0; [ILSMGR. ILS] WD1025. 
DRB0: [ILSMGR. ILS] WD1031. 

RECORD 

5 

STORED 

USING  DRB0: [ILSMGR. ILS] WD1026. 
DRB0: [ILSMGR. ILS] WD1032. 

RECORD 

1 

STORED 

USING  DRB0: [ILSMGR. ILS] WD1027. 
DRB0: [ILSMGR. ILS] WD1032. 

RECORD 

2 

STORED 

USING  DRB0: [ILSMGR. ILS] WD1028. 
DRB0: [ILSMGR. ILS] WD1032. 

RECORD 

3 

STORED 

USING  DRB0: [ILSMGR. ILS] WD1029. 
DRB0: [ILSMGR. ILS] WD1032. 

RECORD 

4 

STORED 

USING  DRB0: [ILSMGR. ILS] WD1030. 
DRB0 : [ ILSMGR. ILS] WD1032 . 

RECORD 

5 

STORED 

$  FIL  1031 

DRB0 : [ ILSMGR. ILS] WD1031 . 

RECORD 

DATA 

12  DK  BLKS,  5  RECORDS 
PRIMARY  FILE 
$  FIL  S1033 

DRB0; IILSMGR.ILSJWD1033.  DOES  NOT  EXIST 

SECONDARY  FILE 
$  OPN  S3 
$  SME  6,2 

DO  YOU  WANT  VARIANCE  RATIOS  PRINTED?  (Y  OR  N) 

->N 

DRB0: [ILSMGR.ILS1WD1033.  RECORD  4  STORED 

$  FIL  S1034 


DRB0 : [ I LSNGR . I LS ] WDl 034.  RECORD  DATA 

12  DK  BLKS,  0  RECORDS 
SECONDARY  FILE 
$  PCO  2,1033 

USER$$DISK: [ILSNGR.ILS]WD1034.  RECORD  5 


USER$$DISK; (ILSMGR.ILS]WD1035.  RECORD  5 
Plot  the  cesaults  using  $  PLR 
$  PLR 

ILS  PLOTTING  ROUTINE 

UP  TO  20  INPUT  FILES 
UP  TO  1024  POINTS  PLOTTED 
ENTER  X-COMPONENT 
->1 

ENTER  Y-COMPONENT 
->2 

ENTER  FILE  NO.,  DISK  NO. 

AND  NO.  CONSECUTIVE  FILES 
DEFAULTS  TO  DRB0 ; [ ILSMGR. ILS] WD1031 . 

->1034, ,2 

SYMBOLS  USED  ARE  A-B 
READ  IN  ANOTHER  FILE?  (Y,N) 

->N 

FILES  BEING  READ 

ENTER  SYMBOL  TO  PLOT,  OR  ENTER  <ALL> 

->ALL 

CC  =  CARTESIAN 
SL  =  SEMI-LOG 
LL  =  LOG- LOG 
->CC 

AUTOMATIC  SCALING?  (Y,N) 

->Y 

GRID?  (Y,N) 

->N 

ENTER  OPTION 
E=EXIT 

N[S,A]=NEW  FILES 

[ST ART, APPEND] 

C=COMPONENTS 
M[I]=MARK  [IDENTIFY] 

R=REASSIGN 

G*GRID 

S» SCALING 

T=TYPE  OF  PLOT 

D=DATA  TO  PLOT 

P[N]=PLOT  [NO  ERASE] 

P[A,E,L]N1,N2= 

FACTOR  ANAL 
H=HARD  COPY 
->E 


STORED 

STORED 


C-5 


